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Abstract 

The  ability  to  fly  satellites  in  close  fonnations  represents  a  capability  that  could 
revolutionize  the  way  satellite  missions  are  designed  in  the  future.  This  study  examines 
three  of  the  primary  formation  flying  designs  and  characterizes  the  affect  that  an 
anomalous  satellite  with  a  slightly  different  cross-sectional  area  would  have  on  the 
stability  of  the  formation.  Following  the  characterization  of  the  effects,  a  controller  is 
implemented  to  mitigate  the  cross-sectional  area  differences  between  the  satellites.  With 
the  addition  of  a  straightforward  controller,  small  cross-sectional  area  differences  can  be 
mitigated  and  corrected  such  that  the  satellites  will  remain  in  close  proximity  and  in  some 
cases  the  fonnation  will  remain  stable. 
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CHARACTERIZING  AND  CONTROLLING  THE  EFFECTS  OF  DIFFERENTIAL 


DRAG  ON  SATELLITE  FORMATIONS 


I.  Background 

The  old  axiom  that  information  is  power  remains  one  of  the  primary  tenets  of 
military  operations.  The  ability  to  observe  ground  events  from  overhead  platforms  can  be 
traced  back  to  the  French  Revolution  when  France  organized  a  company  of  aerostiers,  or 
balloonists  in  the  spring  of  1794  (Richelson,  1999).  Balloon  platforms  were  also  used 
with  limited  success  during  the  American  Revolution,  and  many  scholars  originally 
believed  that  the  airplane’s  primary  wartime  mission  would  revolve  around  aerial 
surveillance  (Richelson,  1999).  When  presented  the  opportunity  to  gather  more 
intelligence  about  their  adversaries  military  leaders  throughout  history  have  jumped  at  the 
idea. 

Upon  the  Soviet  launch  of  Sputnik  on  October  4,  1957,  the  world  of  overhead 
surveillance  changed  forever.  The  Sputnik  launch  led  directly  to  the  creation  of  the 
National  Aeronautics  and  Space  Administration  (NASA)  and  pushed  the  United  States’ 
space  program  into  catch-up  mode.  Fortunately  for  the  United  States,  the  gears  of 
progress  had  already  been  set  in  motion.  On  March  16,  1955,  more  than  two  years  prior 
to  the  launch  of  Sputnik,  the  Air  Force  issued  General  Operational  Requirement  No.  80, 
officially  establishing  a  high-level  requirement  for  an  advanced  reconnaissance  satellite. 
The  document  defined  the  Air  Force  objective  to  be  the  provision  of  continuous 
surveillance  of  "preselected  areas  of  the  earth"  in  order  "to  determine  the  status  of  a 


potential  enemy's  warmaking  capability."  (Richelson,  1999).  The  Air  Force  program, 
originally  Advanced  Reconnaissance  System  (ARS)  then  renamed  SAMOS,  which 
evolved  from  this  directive  experienced  a  number  of  problems  and  delays.  Concern 
about  the  time  schedule  for  SAMOS  prompted  President  Eisenhower  to  approve  the  CIA 
reconnaissance  program  CORONA.  After  12  unsuccessful  launches  the  13th  mission 
finally  returned  a  canister  of  film  from  space  that  was  successfully  recovered  by  the 
United  States  on  August  18,  1960  (Richelson,  1999).  With  the  success  of  CORONA,  the 
high  ground  of  surveillance  had  irreversibly  been  shifted  to  space  borne  platforms. 

Modern  satellite  systems  give  the  United  States  military  an  immense  strategic 
advantage  by  providing  unparalleled  intelligence,  surveillance,  and  reconnaissance  (ISR) 
capabilities.  Traditionally  these  satellites  have  been  large  complex  satellites  that  operate 
independently  or  as  part  of  a  larger  formation  that  provides  near  global  coverage.  These 
complex  systems  require  extra  redundancy  and  extremely  detailed  integration  and  test 
procedures  because  the  failure  of  a  single  subsystem  can  not  only  impact  the  operability 
of  the  individual  satellite,  but  also  endanger  the  mission  of  the  entire  formation. 

However  in  recent  years  mission  planners  have  begun  to  explore  alternative  mission 
profiles  and  satellite  configurations  that  could  reduce  cost,  size,  and  complexity  while 
increasing  imaging  resolution. 

One  idea  that  has  gained  credence  recently  is  the  idea  of  distributing  the 
functionality  of  large  satellites  among  a  group  of  smaller  cooperative  satellites.  This 
concept  involves  flying  a  group  of  satellites  (often  called  a  cluster)  within  very  close 
range  of  each  other  (250m  -  5km)  (Mohammed,  2001:  58).  The  cluster  approach  offers 
mission  planners  much  more  flexibility  in  mission  design  and  lifecycle  planning  because 
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of  the  ability  to  re-position  satellites  in  the  cluster  to  accomplish  different  missions  or  to 
compensate  for  malfunctioning  members  of  the  cluster.  It  offers  additional  redundancy 
and  the  promise  of  graceful  capability  degradation.  Adding  satellites  to  the  cluster  post¬ 
initial  launch  is  also  a  very  real  option  that  allows  for  implementing  new  technologies 
and  enhancing  cluster  performance  (Mohammed,  2001:  58). 

Many  scientific,  military,  and  commercial  space  applications  are  exploring  the  use 
of  clustered  satellites  to  perfonn  distributed  observations  from  space,  synthetic  aperture 
radar  (SAR)  earth  mapping,  magnetosphere  sensing,  interferometry,  and  a  variety  of 
other  missions  (Kitts,  1999:  217).  Many  of  these  missions  are  drawn  to  the  cluster 
approach  because  clusters  can  provide  a  very  large  virtual  aperture.  Since  imaging 
resolution  is  a  function  of  aperture  size,  the  larger  the  aperture  (real  or  virtual)  the  better 
the  resolution.  The  separation  between  instruments  can  be  used  to  create  apertures  that 
allow  for  orders  of  magnitude  improvements  in  ISR  capabilities  (Kitts,  1999:  217). 

In  fact  several  current  missions  are  being  developed  to  validate  cluster  flying 
dynamics,  command  and  control,  and  other  cluster  specific  problems.  The  French  MoD 
recently  launched  a  cluster  of  satellites,  Essaim,  designed  to  demonstrate  the  electro¬ 
magnetic  signal  interception  feasibility  from  space,  and  the  possibilities  of  a  formation 
flying  (swarm)  system  to  prepare  for  coming  fully  operational  systems  (de  Selding, 

2004).  The  Orion  mission  was  developed  to  demonstrate  true  formation  flying  using 
spacebome  carrier  differential  GPS  aboard  low-cost  LEO  satellites.  It  is  scheduled  to  be 
a  single  satellite  flying  in  formation  with  the  two  Emerald  vehicles  (Bauer  et  al,  1999: 
379).  In  addition  many  earth  and  space  science  mission  require  the  use  of  satellite 
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constellations  in  order  to  accomplish  their  mission.  Some  of  these  mission  ideas  will  be 
further  discussed  later  in  the  text. 

There  are  numerous  mission  profiles  that  would  benefit  from  the  use  of  clustered 
satellites.  A  number  of  those  profiles  can  be  accomplished  with  very  loose  formation 
knowledge  (1-100  km).  Unfortunately  on  the  other  end  of  the  spectrum  many  of  the 
possible  science  and  ISR  missions  would  require  much  finer  knowledge  and  control. 
Providing  the  necessary  accuracy  of  micro-meter  and  sometimes  even  pico-meter 
position  knowledge  presents  a  difficult  challenge  for  today’s  command  and  control 
systems  (Baeur  et  al,  1999:  370).  Given  the  large  distances  from  space  down  to  the 
Earth,  it  follows  that  small  errors  in  the  relative  positions  of  the  clustered  satellites  can 
cause  huge  accuracy  and  fidelity  issues  with  the  final  data.  Therefore  one  of  the  primary 
areas  of  research  concerning  clustered  satellites  is  in  the  area  of  command  and  control. 
Currently  there  are  a  number  of  proposed  solutions  to  the  command  and  control  problem, 
but  none  of  them  have  been  successfully  tested  to  the  precision  necessary  to  achieve  pico- 
meter  accuracy.  The  technologies  and  algorithms  for  implementing  these  command  and 
control  schemes  are  still  being  refined,  but  given  the  potential  cost  savings  and 
flexibilities  that  satellite  clusters  could  offer  it  is  only  a  matter  of  time. 

Satellites  in  orbit  experience  a  number  of  perturbations  that  will  test  the  command 
and  control  system  of  any  clustered  formation.  In  low  Earth  orbit  (LEO),  where  most  of 
the  Earth  pointing  systems  will  be  deployed,  one  of  the  primary  perturbations  is 
atmospheric  drag.  Atmospheric  drag  is  a  perturbing  nonconservative  force  that  acts  to 
lower  a  satellite’s  orbit  while  increasing  its  mean  velocity  (Vallado,  2000:  632).  While 
the  acceleration  due  to  drag  does  not  noticeably  effect  satellites  on  an  hour  to  hour  basis, 
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it  is  a  major  perturbation  factor  for  satellites  in  orbits  with  an  altitude  of  less  than  500  km. 
It  becomes  even  more  important  to  understand  when  clustered  satellites  need  to  maintain 
separation  accuracy  on  the  order  of  micro  meters. 


Problem  Statement 

Satellite  clusters  will  ideally  deploy  and  operate  in  an  optimal  manner  every  time. 
However  given  harsh  launch  environments,  excessive  radiation  in  space  and  a  myriad  of 
other  issues  many  satellites  are  not  deployed  in  the  desired  manner  or  develop  problems 
after  a  period  of  operation.  An  anomalous  deployment  or  change  in  the  physical  size  of  a 
satellite  in  a  cluster  can  be  particularly  difficult  to  deal  with  because  of  the  precise 
command  and  control  that  is  necessary  for  the  cluster  to  operate  as  designed.  This  study 
will  investigate  the  effects  of  how  differential  drag  among  satellites  in  a  cluster  formation 
will  affect  the  stability  and  configuration  of  the  cluster. 

Objective 

The  objective  of  this  research  is  to  study  the  effects  of  drag  on  anomalous 
satellites  within  a  cluster  and  determine  how  quickly  the  cluster  fonnation  will 
deteriorate.  This  research  will  examine  the  effects  of  differential  drag  for  satellites  within 
an  in-plane  formation,  an  in-track  formation,  and  a  circular  formation.  In  addition  to 
inspecting  the  effects  of  drag  on  these  formations,  a  controller  will  be  employed  and 
tested  to  detennine  if  a  simple  control  mechanism  may  be  able  to  alleviate  the 
disturbances  created  by  atmospheric  drag  on  the  fonnations. 
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II.  Literature  Review 


Chapter  Overview 

There  has  been  a  surge  of  interest  in  recent  years  regarding  the  enabling 
technologies  for  high  precision  fonnation  flying  satellites.  This  literature  review  will 
examine  three  of  the  major  areas  that  are  distinctive  to  formation  flying  (formation 
design,  formation  navigation  &  control,  and  inter-satellite  communication)  and  the 
current  technology  situation  with  regard  to  those  specific  areas.  In  addition  it  will  briefly 
look  at  the  history  and  current  state  of  modeling  atmospheric  drag. 

Formation  Flying  Overview 

In  addition  to  the  advantages  discussed  in  the  chapter  one,  the  distributed  sensor 
platform  that  formation  flying  promises  is  becoming  more  and  more  enticing  due  to  the 
volume  and  mass  constraints  that  current  launch  vehicle  fairings  place  on  payloads. 

These  constraints  will  restrict  new  monolithic  apertures  to  slightly  larger  or  the  same  size 
as  that  of  the  Hubble  Space  Telescope.  There  are  some  image  resolution  enhancements 
that  can  be  made  by  using  segmented  optics  similar  to  the  optics  planned  for  the  James 
Webb  Space  Telescope,  however  this  approach  is  limited  due  to  structural  dynamics  and 
control  issues  (Leitner,  2004). 
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There  are  three  primary  areas  that  are  distinctive  to  the  formation  flying  problem. 
Those  areas  are  fonnation  design,  formation  navigation  and  control,  and  inter-satellite 
communication  (Leitner,  2004). 

Formation  Design 

Formation  design  can  be  defined  as  the  positioning  of  satellites  into  a  desired 
geometry  to  meet  the  demands  of  the  mission.  It  is  a  difficult  “give  and  take”  problem 
between  optimal  positioning  and  fuel  consumption.  Designing  a  formation  that  maintains 
a  useful  fonnation  with  very  little  or  no  fuel  consumption  is  the  ultimate  goal. 

Sabol  et  al  (2001:  271)  used  the  relative  equations  of  motion  of  two  satellites 
known  as  Hill’s  equations  or  the  Clohessy -Wiltshire  equations  (Chlohessy  and  Wiltshire, 
1960)  to  examine  satellite  formation  design.  They  examined  the  initial  Hill’s  equations 
and  used  a  detailed  derivation  from  Vallado’s  text  (Vallado,  2001 :  374-399)  that  describe 
the  unperturbed  motion  of  two  bodies  (in  close  proximity)  in  circular  orbits  as  a 
preliminary  design  tool. 


2ny-3n2x  =  0 

(2.1) 

y  +  Ink  =  0 

(2.2) 

z  +  n2z  -  0 

(2.3) 

After  solving  analytically  for  these  and  complementing  the  solution  with  numerical 
simulation  of  the  fully  nonlinear  dynamics  with  realistic  force  modeling  they  set  the 
secular  term  to  zero  through  the  constraint 

To  =  -2  x0n  (2.4) 
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“Then  through  algebraic  manipulation  they  showed  that  this  constraint 
results  in  a  displaced  orbit  with  the  same  energy  and  thus  the  same 
semimajor  axis,  as  the  reference  orbit  neglecting  higher-order  terms  in 
eccentricity.  When  the  constraint  is  enforced  Hill’s  equations  become 
(Sabol  et  al,  2001:271)” 

x{t )  =  (x0  / n) sin (nt)  +  x0  cos (nt)  (2.5) 

y(t)  =  (2x0  /  n)  cos  (nt)  -  2x0  sin  (nt)  -  2x0  /  n  +  y0  (2.6) 

z(t)  =  (z0  /  n)  sin(nt)  +  z0  cos  (nt)  (2.7) 

The  above  equations  provided  the  basis  for  Sabol’s  formation  flying  design. 

Three  of  the  formation  flying  designs  that  were  considered  in  their  analysis  were 
the  in-plane  formation,  the  in-track  formation  and  the  circular  formation.  Sabol  et  al 
(2001:  271)  describes  the  in-plane  cluster  design  as  the  simplest  of  all  cluster  designs. 

The  formation  consists  of  a  group  (2+)  of  satellites  occupying  the  same  orbital  plane,  but 
separated  by  mean  anomaly.  This  in-plane  formation  offers  the  advantages  of  a  simple 
design,  straightforward  deployment,  and  uncomplicated  control.  However  this  type  of 
formation  is  limited  in  the  style  of  configurations  it  can  support  and  the  number  of 
satellites  that  can  be  deployed  into  a  reasonable  cluster  (Sabol  et  al,  2001:  271). 
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Figure  1 .  In-plane  formation 

The  in-track  formation  consists  of  satellites  that  all  share  the  same  ground  track. 
At  first  this  concept  seems  identical  to  the  in-plane  design,  and  it  would  be  if  the  Earth 
did  not  rotate.  However  since  the  Earth  does  rotate  satellites  can  only  follow  the  same 
ground  track  if  their  orbital  planes  are  separated  by  a  delta  in  the  right  ascension  of  the 
ascending  node,  which  compensates  for  the  rotating  Earth.  This  formation  offers  the 
advantage  of  identical  ground  tracks,  but  will  also  have  a  difficult  time  supporting 
multiple  configuration  shapes  with  numerous  satellites  (Sabol  et  al,  2001:  272). 


i» 

Figure  2.  In-track  formation  (Sabol  et  al,  2001 :  272) 
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The  final  type  of  formation  that  this  study  examines  is  the  circular  formation.  In 
the  circular  formation  the  satellites  maintain  a  constant  distance  from  each  other.  For  two 
arbitrary  points  in  the  cluster  their  will  most  likely  be  differences  in  inclination,  right 
ascension  of  the  ascending  node,  argument  of  perigee,  and  mean  anomaly.  Sabol  et  al 
(2001 :  272)  give  the  following  insight  into  the  circular  fonnation  design  in  terms  of 
Keplerian  mechanics, 

“Consider  a  circular  reference  orbit  inclined  at  90  degrees  with  a 
satellite  at  the  equator  (ascending  node).  Now  consider  another 
satellite  in  a  slightly  elliptical  orbit  also  on  the  equator  and 
separated  from  the  reference  orbit  by  a  small  amount  of  ascending 
node.  The  second  satellite  is  at  apogee  and,  therefore  will  fall 
behind  the  reference  satellite  as  they  both  proceed  towards  the 
North  Pole.  Because  both  orbits  are  polar,  the  satellites’  paths  will 
cross  at  the  poles,  but  the  reference  satellite  reaches  the  pole  first. 

On  the  other  side  of  the  North  Pole,  the  second  satellite  is  lower  in 
altitude  than  the  reference  satellite  and  begins  to  catch  up.  Both 
satellites  reach  their  descending  nodes  at  the  same  time,  and  the 
second  satellite  is  now  at  perigee.  The  second  satellite  continues  to 
advance  ahead  of  the  reference  satellite  and  reaches  the  South  Pole 
first,  where  the  paths  cross  once  again.  On  the  other  side  of  the 
South  Pole,  the  reference  satellite  begins  to  catch  up  as  the  second 
satellite  increases  in  altitude  towards  its  apogee  and  ascending 
node.” 
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The  circular  formation  offers  two  characteristics  make  it  attractive  to  cluster  designers. 
One  advantage  is  that  the  two  satellites  maintain  a  constant  distance  from  each  other  at  all 
times.  The  other  advantage  that  is  unique  to  the  circular  fonnation  is  that  unlike  the  in¬ 
plane  and  in-track  formations  the  circular  fonnation  presents  a  two-dimensional  array 
(Sabol  et  al,  2001:  272-273). 

Formation  Navigation  and  Control 

The  implementation  of  formation  flying  satellites  acting  in  hannony  to  achieve 
military  and  scientific  Earth  sensing  objectives  will  require  stringent  control  thresholds. 
Formation  control  is  responsible  for  rejecting  disturbances  and  keeping  the  formation  in 
the  desired  geometry.  Formation  maintenance  and  control  at  the  necessary  level  is  one  of 
the  more  difficult  problems  facing  the  widespread  deployment  of  clustered  satellites. 
Formation  control  and  relative  navigation  pose  some  very  interesting  challenges  in  the 
areas  of  (Bauer  el  al,  1999:  376): 

1 .  Onboard  sensing  of  relative  and  absolute  vehicle  position/attitudes. 

2.  Maneuvering,  retargeting,  collision  avoidance,  and  aperture  optimization 
including  resource/task  allocation  within  the  fleet. 

3.  Modeling  the  orbital  mechanics  and  the  impact  of  differential  drag  and  solar 
disturbances 

4.  Fleet  and  vehicle  autonomy,  including  high-level  fault  detection  and  recovery  to 
enhance  mission  robustness 
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5.  Decentralized  control  and  computation  for  a  fleet  of  many  (e.g.  from  16  to 
hundreds)  vehicles. 

6.  Testbeds  and  simulations  to  validate  the  various  sensing  and  control  concepts. 

Many  current  spacecraft  systems  rely  on  ground  based  command  and  control  to 
ensure  proper  positioning  and  attitude  of  individual  spacecrafts.  However  nearly  all  of 
those  systems  are  single  satellites  that  are  widely  spaced  and  don’t  depend  on  the  relative 
positioning  data  between  spacecraft.  When  dealing  with  a  clustered  group  of  satellites 
ground  based  command  and  control  becomes  much  more  cumbersome.  Attitude  and 
position  data  are  needed  instantaneously.  Ground-based  systems  would  be  very  intricate, 
weighed  down  by  data,  and  possibly  unable  to  provide  timely  corrective  commands. 
Thus,  a  number  of  institutions  have  been  focusing  on  developing  autonomous  satellite- 
based  command  and  control  algorithms  and  schema  to  facilitate  clustered  satellite  flight 
(Bauer  et  al,  1999:  376). 

Ideally,  autonomous  formation  flying  involves  taking  continuous  position  and 
velocity  data  that  determine  the  state  of  the  formation.  Those  measurements  are  then 
used  to  maintain  the  current  configuration  or  to  transition  to  a  new  formation  style 
without  the  application  of  any  outside  controls.  The  configuration  of  the  cluster  includes 
not  only  the  distances  between  all  pairs  of  spacecraft  in  the  cluster,  but  also  the 
orientation  of  the  satellites  in  a  coordinate  frame  defined  by  the  array’s  internal  geometry 
(Bauer  et  al,  1999:  376).  Initializing,  targeting,  and  maneuvering  will  all  task  the 
autonomous  command  and  control  system  in  different  manners  such  that  the  control 
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system  must  be  designed  and  implemented  to  switch  between  the  various  system  models 
and  control  schemas. 

Cooperative  formation  control  can  be  achieved  in  several  manners.  A  typical 
control  design  for  a  small  number  of  satellites  (16  or  less)  in  formation  is  to  employ  a 
master/slave  control  structure.  Master/slave  scenarios  involve  a  single  spacecraft  that 
acts  as  the  leader  (master)  and  issues  commands  to  the  other  spacecrafts  (deputies  or 
slaves).  This  can  be  accomplished  by  actual  commands  being  sent  to  the  other  members 
of  the  cluster  or  by  the  slaves  autonomously  reacting  to  the  maneuvers  of  the  master  in 
order  to  maintain  the  proper  configuration  (Bauer  et  al,  1999:  376).  This  master/deputy 
control  structure  works  well  with  a  small  number  of  satellites,  and  can  be  scaled  up 
slightly  by  dividing  the  cluster  into  a  number  of  subsets.  However  as  the  number  of 
satellites  increase  the  overhead  becomes  unbearable  and  the  system  begins  to  break 
down,  necessitating  a  distributed  control  system  when  the  formation  is  large  (Bauer  et  al, 
1999:  376). 

A  decentralized  control  system  demands  that  each  node  in  the  system  process  its 
own  measurement  data  in  parallel  with  all  of  the  other  nodes.  Detected  failures  then  can 
be  mitigated  so  that  the  cluster  performance  degrades  gracefully  (Bauer  et  al,  1999:  376). 
However  none  of  these  navigation  and  control  systems  have  been  tested  and  verified  in  a 
relative  environment  or  for  the  desired  amount  of  time  that  fonnations  would  stay  on 
orbit. 
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Inter-Satellite  Communication 


Along  with  precise  navigation  and  control  requirements,  the  communication  link 
between  the  formation  flyers  is  critical  to  success.  The  data  buses  must  be  extremely 
robust  and  reliable.  There  remains  a  substantial  amount  of  work  to  be  done  in  defining 
requirements  for  communication  bandwidth,  time  synchronization,  and  the  transfer  of 
precision  control  commands  (Francis  et  al,  2003). 

The  importance  of  this  system  is  akin  to  the  mailman  in  the  postal  system. 
Everything  else  may  go  perfectly  but  nobody  receives  the  mail  without  the  mailman.  The 
formation  control  laws  will  all  be  implemented  through  the  communication  system,  so  a 
lack  of  integrity  in  the  system  will  mean  a  loss  of  control  that  could  be  catastrophic  even 
if  it  happens  for  only  a  few  seconds.  This  is  especially  true  if  the  formation  is  to  employ 
a  decentralized  command  structure  that  relies  heavily  on  each  communication  unit 
functioning  flawlessly. 

State  of  the  Technology 

Formation  Design 

Sabol  et  al  (2001:  274)  also  examined  the  relative  stability  of  two  identical 
satellites  initially  separated  by  1  km,  placed  into  in-plane  fonnations  and  in-track 
formations.  They  looked  at  two  different  cases  for  each  fonnation.  The  initial  case  was  a 
100-revolution-per-seven-nodal-day  repeat  ground  track  cycle  and  the  second  situation 
involved  a  14  revolution  per  nodal  day  repeat  ground  track  cycle.  They  then  propagated 
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the  orbit  scenarios  using  the  DSST  Averaged  Orbit  Generator  (AOG)  which  includes  the 
following  modeled  perturbations:  geopotential,  atmospheric  drag,  luni-solar  third  body 
point  mass  effects,  and  solar  radiation  pressure  (Sabol  et  al,  2001 :  274). 

The  in-plane  formation  proved  to  be  extremely  stable  in  the  (100:7)  case  and  only 
showed  a  diversion  of  a  few  meters  over  a  year  long  timeframe.  During  the  daily  repeat 
ground  track  (14: 1)  tesseral  resonance  had  a  much  greater  effect  and  caused  the  satellites 
to  drift  nearly  0.2  km  closer  together  over  a  year.  Over  time  (the  next  year)  the  satellites 
should  drift  back  to  near  their  original  positions  (Sabol  et  al,  2001 :  274). 

The  in-track  formation  satellites  share  the  same  ground  tracks  so  they  encounter 
the  same  gravitational  effects.  However  the  two  satellites  are  in  slightly  different  orbital 
planes  so  the  high-density  atmosphere  simulation  causes  a  slight  along-track  drift 
between  the  two  satellites.  This  drift  can  be  shown  to  be  on  the  order  of  about  50m  after 
one  year  (Sabol  et  al,  2001:  274). 

Both  of  these  simple  formation  styles  offer  promising  stable  designs  that  should 
require  minimal  station  keeping  propulsion  and  maintenance.  However  these  studies  all 
examined  perfectly  identical  satellites  in  the  fonnations,  any  slight  differences  in  the 
satellites  would  produce  different  results. 

The  circular  formation  was  also  studied,  albeit  in  a  slightly  different  manner 
(Sabol  et  al,  2001 :  274).  Three  satellite  trajectories  were  generated  for  the  circular 
cluster.  Satellite  one  represents  the  center  of  the  circular  fonnation  and  the  reference 
satellite.  Satellites  two  and  three  are  on  the  circle  phased  at  270  degrees  and  180  degrees. 
This  design  puts  satellite  two  directly  in  front  of  satellite  one  in  the  in-track  direction  and 
satellite  three  to  the  right  in  the  negative  normal  (cross-track)  direction.  For  their  study 
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Sabol  et  al  (2001 :  274)  used  an  orbit  of  800  km  and  began  propagating  the  orbit  scenarios 
using  the  DSST  AOG.  They  found  the  propagations  in  the  presences  of  perturbations 
showed  that  the  circular  formation  was  unstable  at  an  altitude  of  800  km.  The  primary 
disrupting  factor  was  found  to  be  the  Earth’s  oblateness  or  the  J2  effect.  The  J2  effect 
disturbs  the  formation  in  two  primary  manners:  1)  through  the  differential  precession  of 
the  orbital  planes  and  2)  through  the  shifting  of  the  argument  of  perigee  (Sabol  et  al, 
2001:274-275). 

Because  satellites  one  and  two  have  slightly  different  inclinations  the  secular  J2 
effect  causes  the  right  ascension  of  the  ascending  node  for  the  different  orbits  to  precess 
at  different  rates.  This  results  in  the  orbital  planes  drifting  apart  and  error  growth  in  the 
nonnal  (or  cross-track)  direction.  This  error  growth  expands  fairly  quickly  and  after  four 
days  in  the  800  km  orbit  the  separation  between  satellite  one  and  two  could  increase  by 
up  to  25%  (Sabol  et  al,  2001:  276).  For  the  worst  case  circular  fonnation,  formation 
maintenance  can  cost  up  to  38  m/s  per  year  per  kilometer  of  formation  radius  (Sabol  et  al, 
2001:  276). 

Rotation  of  the  argument  of  perigee  is  another  major  disruption  on  the  circular 
formation.  This  effect  is  independent  of  formation  size,  but  the  budget  for  maintaining 
constant  argument  of  perigees  is  nearly  1 1  m/s  per  year  (Sabol  et  al,  2001 :  276). 

Given  the  major  disrupting  force  that  the  J2  effect  causes  upon  circular  formations 
it  will  be  necessary  to  make  fonnation  keeping  maneuvers  on  a  near  daily  basis  (Sabol  et 
al,  2001:276). 
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Formation  Navigation  and  Control 

NASA/GSFC  is  currently  pursuing  a  decentralized  approach  that  involves  using  a 
stand-alone  GPS  point  solution  to  maintain  the  spacecraft  formation  (Carpenter  et  al, 
2003:  3).  In  this  decentralized  approach  each  satellite  transmits  and  receives  data  to  and 
from  each  of  the  other  satellites  in  the  formation.  This  allows  relative  states  to  be 
computed  without  the  need  for  direct  measurement  of  the  inter-satellites  states.  If  enough 
processor  capacity  is  available  the  GPS  measurement  data  could  be  processed  for 
improved  accuracy.  The  addition  of  instruments  that  would  allow  one  or  more  of  the 
vehicles  to  take  relative  measurements  between  itself  and  the  other  formation  members 
would  provide  additional  data  that  could  maximize  the  relative  navigation  accuracy.  This 
particular  approach  with  the  multi-faceted  measurements  will  demand  a  lot  of  processing 
power.  In  addition  an  approach  similar  to  this  may  be  necessary  to  obtain  the  necessary 
pico-meter  accuracy  that  will  be  necessary  for  a  number  of  the  desired  satellite  cluster 
missions  (Carpenter  et  al,  2003:  3). 

Stanford  University  has  also  been  developing  methods  for  centralized  and 
decentralized  control  of  a  satellite  cluster.  Their  team  has  demonstrated  several 
estimation  architectures  that  could  be  used  in  a  differential  carrier-phase  GPS  relative 
sensing  for  larger  (>16)  cluster  of  satellites  (Corazzini  et  al,  1997).  They  have  used  the 
Formation  Flying  Testbed  at  Stanford  to  demonstrate  a  multi-level  cluster  control  system 
that  includes  a  coordinator,  a  planner,  and  distributed  regulators  using  a  high  fidelity  orbit 
simulator.  A  linear  programming  approach  allows  their  system  to  rapidly  solve  for  the 
optimal  formation  maneuvers  using  linearized  group  dynamics.  A  portion  of  their 
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experimentally  analyzed  autonomous  formation  flying  algorithms  will  be  demonstrated 
on  Orion  (How  et  al,  1998). 

Currently  the  NASA  Goddard  Space  Flight  Center  Earth  Observing  1  (EO-1)  is  on 
orbit  and  employing  a  manner  of  fonnation  flying.  The  EO-1  formation  flying 
demonstration  is  controlled  from  the  ground  (not  cross-linked)  using  magnetospheric 
measurements.  However  it  is  only  important  that  the  satellites  be  arranged  in  a  particular 
shape  (tetrahedron)  at  a  single  point  (apogee)  in  the  orbit  (Folta  et  al,  2002).  Similarly 
the  Solar  Imaging  Radio  Array  (SIRA)  mission  will  require  only  a  loose  control  of 
positions  and  spacing  such  that  all  spacecraft  are  within  a  spherical  region  (MacDowell  et 
al  2005).  These  two  missions  are  both  formation  flying  systems,  but  their  controlling 
schemes  are  loose  enough  that  controls  need  to  be  applied  only  once  or  twice  an  orbit. 

In  addition,  a  number  of  other  agencies  are  developing  fonnation  flying 
technologies  and  base  lining  clusters  for  use  in  future  missions.  The  European  Space 
Agency  has  developed  a  sophisticated  cluster  fonnation  to  study  the  Earth’s 
magnetosphere  (Roux,  1998).  The  low-cost  technology  demonstration  Orion-Emerald 
program  developed  by  a  group  of  universities  will  demonstrate  the  use  of  carrier-phase 
differential  GPS  as  a  primary  sensor  for  formation  flying  (Kitts  et  al,  1999)(How  et  al, 
1998). 

Air  Force  Research  Laboratory’s  (AFRL)  technology  demonstration  missions 
TechSat  21  attempted  to  tackle  some  of  the  extreme  system  level  challenges  and 
technology  hurdles  that  face  future  formation  flying  missions.  Their  efforts  uncovered 
even  more  technology  challenges  and  eventually  led  to  the  program  cancellation  (Cobb, 
2005). 
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For  many  of  the  envisioned  formation  flying  missions  continuous  six  degrees  of 
freedom  formation  control  must  be  implemented  at  levels  that  are  orders  of  magnitude 
more  precise  than  any  spacecraft  on  orbit  today.  “This  will  necessitate  propulsion 
systems  with  constantly  varying  levels  of  thrust,  unlike  impulsive  systems  currently  in 
use  today,”  (Leitner,  2004). 

Inter-Satellite  Communication 

Current  satellite  communication  systems  are  primarily  designed  for  uplink  and 
downlink  purposes.  Inter-satellite  communication  offers  a  different  type  of  challenge  that 
will  require  increased  precision  and  robustness.  A  full  analysis  of  the  effect  of  variable 
distance  communication  delays  between  formation  flying  satellites  and  how  it  will  affect 
closed-loop  control  performance  must  be  completed  and  evaluated  (Leitner,  2004). 

Modeling  Atmospheric  Drag  (Vallado,  2001:  524-537) 

Satellites  orbit  around  the  Earth  based  upon  the  laws  of  two-body  motion 
described  by  Kepler.  During  their  orbits  they  experience  a  number  of  perturbations  that 
affect  the  basic  equations  of  two-body  motion;  four  of  the  largest  perturbations  are 
atmospheric  drag,  third-body  effects,  Earth  oblateness,  and  solar  radiation  pressure.  The 
relative  importance  of  these  effects  varies  greatly  based  upon  the  height  of  the  satellite’s 
orbit.  In  low  Earth  orbit  (LEO)  up  to  an  altitude  of  about  500  km,  atmospheric  drag  is  the 
primary  perturbation  force.  The  Earth’s  oblateness  also  causes  a  pronounced  departure 
from  two  body  motion  and  is  the  second  most  pronounced  perturbation  for  LEO 


19 


satellites.  When  satellites  move  into  orbits  above  the  LEO  range,  solar  radiation  pressure 
and  third-body  effects  begin  to  play  a  dominant  effect  on  the  orbital  perturbations. 

Atmospheric  drag  is  caused  by  atmospheric  particles  colliding  with  the  satellite 
and  impeding  the  satellite’s  motion.  The  acceleration  due  to  the  drag  force  changes 
based  upon  the  density  of  the  particular  atmosphere,  and  it  can  be  expressed  as: 


2  m 


(2.8) 


The  coefficient  of  drag,  Cd ,  is  a  dimensionless  quantity  that  expresses  the  satellite’s 

susceptibility  to  drag  forces.  The  exposed  cross-sectional  area,  A,  is  the  area  of  the 
satellite  that  is  normal  to  the  velocity  vector  of  the  satellite.  It  can  be  a  difficult  quantity 
to  estimate  accurately  unless  the  configuration  and  attitude  of  the  satellite  is  precisely 
known.  The  atmospheric  density,  p  ,  comes  from  one  of  the  numerous  atmosphere 
models  that  are  available,  but  is  difficult  to  predict  because  of  it’s  variable  nature.  All 
three  of  these  quantities  are  difficult  to  accurately  measure  and  predict  so  the  science  of 
determining  precise  accelerations  due  to  atmospheric  drag  still  suffers  from  some 
ambiguity. 

To  accurately  predict  the  effects  of  drag  on  satellites  in  specific  orbits  we  need  to 
be  able  to  accurately  model  the  atmosphere  at  the  altitude  of  interest.  Unfortunately  the 
Sun’s  interaction  with  the  upper  atmosphere  and  the  Earth’s  magnetic  field  cause  the 
properties  of  the  atmosphere  to  shift  and  change  such  that  there  will  always  be  a  level  of 
uncertainty  while  modeling  the  Earth’s  atmosphere. 
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Despite  the  known  difficulties,  scientists  and  astronomers  have  developed 
numerous  atmospheric  density  models  over  the  past  twenty  years  in  order  to  meet  the 
changing  accuracy  demands  of  specific  projects.  Most  of  these  models  include  some 
manner  of  exponential  decay  as  altitude  goes  up.  Those  models  can  be  separated  into  two 
main  groups:  static  models  and  time-varying  models.  Gaposhkin  and  Coster  (1998) 
discuss  many  of  the  atmospheric  density  models  in  detail  and  come  to  the  conclusion  that 
no  model  is  best  for  all  applications.  Computing  power  and  accuracy  requirements  are 
the  primary  deciding  factors  when  choosing  the  best  drag  model  for  a  particular 
application.  However  more  computing  power  does  not  necessarily  mean  better  accuracy 
and  other  factors  such  as  length  of  simulation  and  sections  of  atmosphere  pertinent  to  the 
given  experiment  need  to  be  considered  when  determining  which  atmospheric  model  best 
fits  the  application. 
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III.  Methodology 


Chapter  Overview 

The  objective  of  this  research  is  to  study  the  effects  of  drag  on  anomalous 
satellites  within  a  cluster  and  determine  how  quickly  the  cluster  fonnation  will 
deteriorate.  This  research  will  examine  the  effects  of  differential  drag  for  satellites  within 
an  in-plane  formation,  an  in-track  formation,  and  a  circular  formation.  In  addition  to 
inspecting  the  effects  of  drag  on  these  formations,  a  controller  will  be  employed  and 
tested  to  detennine  if  a  simple  control  mechanism  may  be  able  to  alleviate  the 
disturbances  created  by  atmospheric  drag  on  the  fonnations. 

The  selected  approach  for  this  research  is  to  simulate  all  of  the  cluster  formations 
with  two  satellites  in  the  formation.  This  approach  reduces  the  complexities  of  dealing 
with  numerous  satellites,  but  allows  detailed  study  of  the  anomalous  satellite  and  its 
position  relative  to  a  control  satellite.  One  of  the  satellites  in  the  formation  is  considered 
the  control  satellite  and  its  orbital  characteristics  will  remain  constant  throughout  all  of 
the  experiments  (unless  explicitly  stated).  The  second  satellite’s  (the  deputy’s)  orbital 
elements  will  change  based  upon  the  specific  scenario  that  is  being  examined. 

Control  Satellite 

Throughout  the  different  test  cases,  the  control  satellite  maintains  constant  initial 
orbital  elements  (excluding  true  anomaly).  It  is  orbiting  the  Earth  in  a  circular  orbit 
(eccentricity  =  0)  at  an  altitude  of  250  kilometers  in  a  polar  orbit  (inclination  =  x  12 
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radians).  The  control  satellite  has  a  right  ascension  of  the  ascending  node  of  zero  radians 
and  the  simulation  starts  when  the  control  satellite  is  on  the  equator  (argument  of  latitude 
=  0). 


Deputy  Satellite 

The  deputy  satellite’s  orbital  elements  will  change  for  each  cluster  formation 
scenario  that  is  run.  In  addition  to  the  changes  in  orbital  elements  that  detennine  the 
specific  formation  style,  the  deputy  satellite  will  also  experience  anomalies  that  will 
cause  it  to  have  a  different  cross  sectional  area  (and  thus  experience  a  different 
acceleration  due  to  drag)  than  the  control  satellite. 


Atmospheric  Drag  Modeling 

Satellites  in  low  orbits  are  constantly  fighting  against  the  forces  of  drag,  and  the 
lower  the  orbit  altitude  the  stronger  the  forces  of  drag.  Because  satellites  are  fighting 
against  the  air  particles  to  move  forward,  atmospheric  drag  decreases  the  energy  of  the 
satellite’s  orbit.  Therefore,  the  acceleration  due  to  drag  is  a  negative  (opposite  of  the 
satellite’s  velocity  vector)  acceleration  and  it  can  be  expressed  as: 


Ad 


1  (CaA)  2 

- - Pv 

2  m 


(3.1) 


Where: 

•  Ad  is  the  acceleration  due  to  drag 
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Cd  is  the  coefficient  of  drag 


• 

•  ^4  is  the  cross  sectional  area  perpendicular  to  velocity  vector 

•  m  is  the  mass  of  the  satellite 

•  p  is  the  density  of  the  atmosphere  at  position  of  the  satellite 

•  v  is  the  magnitude  of  the  velocity 

The  cross  sectional  area  perpendicular  to  the  velocity  vector  can  be  estimated  if  you 
know  the  attitude  of  the  satellite  in  relation  to  the  direction  of  travel.  For  this  study  a 
cross  sectional  area  of  2  square  meters  was  used  for  the  control  satellite.  This  simulates  a 
small  to  medium  sized  satellite  with  solar  arrays  deployed.  This  study  also  used  a  cross 
sectional  area  of  1.5  meters  squared  to  simulate  a  satellite  experiencing  an  anomaly  such 
as  a  partially  deployed  solar  array.  At  times  in  this  research,  a  much  larger  cross 
sectional  area  difference  will  be  used  to  demonstrate  the  effects  of  drag  in  a  timelier 
manner.  The  coefficient  of  drag  is  usually  determined  experimentally  and  is  based  upon 
the  complex  dependencies  of  shape,  attitude,  flow  conditions,  and  spacecraft  drag.  This 
study  uses  a  number  of  2.2,  which  is  a  typical  value  for  satellites  (Larson  and  Wertz, 

1991 :  143).  This  research  uses  an  arbitrary  mass  of  100  kilograms  for  each  satellite  in 
the  simulation.  This  mass  is  the  upper  limit  for  micro-satellites. 

As  discussed  in  the  literature  review  there  are  numerous  different  methods  for 
modeling  the  atmosphere  and  computing  atmospheric  density.  For  this  study  a  relatively 
simple,  but  accurate,  exponential  atmospheric  model  will  be  used  (Vallado,  2001:534- 
535).  This  simple  model  assumes  the  density  of  the  atmosphere  decays  exponentially 
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with  increasing  altitude.  It  also  assumes  a  spherically  symmetrical  distribution  of 
particles,  in  which  the  density,/?  ,  varies  exponentially  according  to  Vallado  (2001:  535): 


p  =  paEXP[-—, 7^4 

ti 


(3.2) 


Where: 

•  p0  is  the  reference  density 

•  h0  is  the  reference  altitude 

•  hellp  is  the  actual  altitude  above  the  Earth 

•  H  is  the  scale  height 

By  using  the  following  table  to  determine  values  for  specific  bands  of  the  atmosphere,  the 
exponential  atmospheric  model  can  become  quite  accurate  in  predicting  the  nominal 
density  at  a  given  altitude. 


25 


Altitude 

hellp 

(km) 

Base 
Altitude 
h0  (km) 

Nominal 

Density 

p„(kg/m3) 

Scale 
Height 
H  (km) 

Altitude 

hellp 

(km) 

Base 
Altitude 
ha  (km) 

Nominal 

Density 

pa  (kg/m3) 

Scale 

Height 

H  (km) 

0-25 

0 

1.225 

7.249 

150-180 

150 

2.070  x  IO"9 

22.523 

25-30 

25 

3  899  x 10‘2 

6.349 

180-200 

180 

5.464  x  IO'10 

29.740 

30-40 

30 

1.774  x  IO'2 

6.682 

200-250 

200 

2.789  x  IO'1" 

37.105 

40-50 

40 

3 972 x 10"3 

7.554 

250-300 

250 

7.248  x  10" 11 

45.546 

50-60 

50 

1.057  x  10-3 

8  382 

300-350 

300 

2.418  x  10~“ 

53.628 

60-70 

60 

3.206  x  !0~* 

7.714 

350-400 

350 

9.518  x  IO'12 

53.298 

70-80 

70 

8.770  x  IO'5 

6.549 

400-450 

400 

3.725  x  10"12 

58.515 

80-90 

80 

1.905  x  10' 5 

5.799 

450-500 

450 

1.585  x  IO-12 

60.828 

90-100 

90 

3.396  x  IO-6 

5.382 

500-600 

500 

6  967  x  IO-13 

63  822 

100-110 

100 

5.297  x  I0‘7 

5.877 

600-700 

600 

1.454  x  10" 13 

71.835 

110-120 

no 

9  661  x  10"8 

7.263 

700-800 

700 

3.614  x  IO-14 

88.667 

120-130 

120 

2.438  x  lO-8 

9.473 

800-900 

800 

1.170x  IO-14 

124.64 

130-140 

130 

8.484  x  10'9 

12636 

900-1000 

900 

5.245  x  IO"15 

181.05 

140-150 

140 

3.845  x  10~9 

16  149 

1000- 

1000 

3.019  x  IO'15 

268.00 

Figure  3.  Exponential  atmosphere  modeling  table  (Vallado,  2001:  537) 

This  atmospheric  density  modeling  approach  works  sufficiently  well  for  design  level 
studies  such  as  the  one  being  conducted;  however,  highly  accurate  studies  might  choose 
to  use  more  sophisticated  and  accurate  models.  The  atmospheric  density  derived  from 
the  above  equations  and  tables  will  then  be  incorporated  into  equation  3.1  and  assimilated 
into  the  orbit  propagator. 

Orbit  propagation 

This  study  uses  an  orbit  propagator  based  solely  on  two  body  motion  and  the 
effects  of  drag.  Third  body  and  non-spherical  Earth  perturbations  (J2  effect)  were  not 
taken  into  account  because  the  focus  of  this  work  is  to  understand  how  differential  drag 
will  affect  satellite  formations.  The  primary  effects  these  perturbations  have  upon  orbits 


26 


are  shifting  the  longitude  of  the  ascending  node  and  moving  the  argument  of  perigee 
(Wiesel  1997:  88).  While  important,  it  will  not  significantly  effect  the  relative  position  of 
satellites  in  an  in-plane  or  an  in-track  cluster  (Sabol  et  al,  2001 :  274-275).  The  J2  effect 
changes  the  ascending  node  ( Q  )  by  causing  it  to  precess  at  the  rate  (Wiesel,  200 1:88) 


Q  = 


3nJ2RE2 
2a2 (\  -e2  )1 


cos  i 


(3.3) 


where: 

•  n  is  the  mean  anomaly 

•  J2  is  a  dimensionless  number  that  characterizes  the  departure  of  a  body  from  a 
true  sphere,  for  the  Earth  J2  =0.001082 

•  Re  is  the  radius  of  the  Earth 

•  a  is  the  satellite’s  semi-major  axis 

•  e  is  the  eccentricity  of  the  orbit 

•  i  is  the  inclination  of  the  orbit 

In  addition  the  J2  effect  causes  elliptical  orbits  to  rotate  in  their  own  plane  at  a  rate  of 
(Wiesel,  2001:  88): 


d>  =  - 


3nJ2R2 

la2(\-e2)2 


(—sin^  i-2) 
2 


(3.4) 
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As  previously  discussed  in  the  literature  review  and  to  be  discussed  further  later  in  this 
chapter  both  the  in-plane  and  in-track  formations  offer  formations  where  all  of  the  orbits 
of  the  satellites  in  the  cluster  have  the  same  inclination,  semi-major  axis,  and  eccentricity. 
This  will  cause  Q  and  6)  to  remain  equal  for  the  satellites  clustered  in  and  in-plane  and 
in- track  formations. 

However  the  J2  effect  does  play  a  significant  role  in  the  stability  of  satellites  in  a 
circular  cluster  (Sabol  et  al,  2001:  274).  Equations  3.3  and  3.4  show  that  satellite  orbits 
with  different  inclinations  and  eccentricities  will  experience  different  Q  and  co .  As 
discussed  in  the  literature  review  the  J2  effect  is  the  primary  reason  that  the  circular 
formation  is  unstable  and  a  two  satellite  circular  formation  may  require  up  to  50  m/s  per 
year  per  kilometer  of  separation  and  daily  orbital  corrections  to  maintain  a  stable 
formation  (Sabol  et  al,  2001 :  276).  This  particular  study  is  done  under  the  context  that  a 
controller  of  some  type  is  already  in  place  that  mitigates  the  effect  of  the  non-spherical 
Earth  (and  other  less  disruptive  perturbations)  on  the  circular  formation.  Without  the 
controller  that  accomplishes  this  task  the  circular  formation  is  unstable  from  the  very 
beginning  and  cannot  be  considered  a  feasible  long  tenn  formation.  Thus  it  would  not  be 
important  how  differential  drag  would  effect  the  formation.  In  addition,  solar  radiation 
perturbations  were  not  considered  because  of  the  low  relative  magnitude  of  the 
perturbation  in  comparison  with  atmospheric  drag  (Vallado,  2001 :  646).  The  differential 
equations  expressed  in  the  Earth  centered  inertial  (ECI)  frame  used  in  the  propagator  are 
as  follows: 


x  =  -(///r3)x 


2  m 


(3.5) 
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■■  ,  ,  1  (C,A)  , 

y=  (ju/r)y  v  pv  y/v 

(3.6) 

2  m 

..  ,  .  2.  l(M  2  -  i 

z  =  -{plr  )z - - - pv  z/v 

2  m 

(3.7) 

For  the  initial  cases  that  examine  the  effects  of  drag  on  the  formation,  the  orbits  of  the 
control  satellite  and  the  deputy  satellite  are  propagated  independently  and  then  compared 
at  specific  time  intervals  (every  10  seconds)  to  determine  the  difference  in  the  orbits. 

Initial  Conditions  Conversion 

The  MATLAB  program  written  for  this  research  accepts  the  classical  orbital 
elements  as  input  parameters.  However  in  order  to  propagate  the  orbits  of  the  satellites 
we  need  to  convert  from  the  initial  classical  orbital  elements  to  an  ECI  frame.  This  is 
done  by  first  converting  the  traditional  orbital  elements  into  the  perifocal  coordinate 
system  and  then  completing  a  coordinate  transformation  that  transforms  them  into  the 
ECI  frame. 

Convert  from  Orbital  Elements  to  ECI  Frame 

For  the  chosen  equations  of  motion  the  coordinates  of  the  satellites  must  be 
expressed  in  x,  y,  and  z  components  in  the  ECI  frame.  The  MATLAB  program  written 
for  this  research  accepts  the  traditional  six  orbital  elements;  inclination  (i),  right 
ascension  of  the  ascending  node  (Q),  eccentricity  (e),  argument  of  perigee  ( co ),  semi¬ 
major  axis  (a),  and  either  mean  anomaly  or  true  anomaly  ( v)  and  converts  these  elements 
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to  position  and  velocity  vectors  in  the  perifocal  coordinate  system.  The  perifocal  position 
and  velocity  vectors  are  obtained  by  the  following  equations  (Wiesel,  1997:  65) 


perifocal 


cos  u/(l  +  e  cosy) 
sin  u/(l  +  e  cost/) 
0 


(3.8) 


v 


perifocal 


-y/plp 
-jp!  Pie) 
0 


sinu 
+  cos  v 


(3.9) 


where: 

•  p  is  the  orbit  semi-parameter,  =  a(\  -e2) 

•  p  is  the  gravitational  parameter  of  the  Earth 

Then,  those  perifocal  position  and  velocity  vectors  go  through  a  coordinate 
transformation  that  transforms  them  from  perifocal  coordinates  to  ECI  coordinates.  That 
coordinate  transformation  can  be  expressed  as: 


v  —  /?  R  R  r 

eci  VQ 1  S' 1  at r  perifocal 

(3.10) 

^ eci  R'Q.^i^co^ perifocal 

(3.11) 

Where  the  rotation  matrices  are  as  follows: 
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(3.12) 


R 


0) 


cos(-<x>)  sin(-ft))  0 
-sin(-<z>)  cos(-<z>)  0 
0  0  1 


R.  = 


0  0 
cos(-z')  sin(-z) 
-sin(-z')  cos  (-/) 


(3.13) 


R 


n 


cos(-Q)  sin(-Q)  0 
-sin(-Q)  cos(Q)  0 
0  0  1 


(3.14) 


ECI  to  RTN 

Traditionally  two  satellite  motion  has  been  expressed  in  the  Radial,  Tangential, 
Normal  (RTN)  frame  with  one  (the  control  satellite  in  this  case)  of  the  satellites 
representing  the  origin  of  the  coordinate  system.  In  order  to  transform  coordinates 
expressed  in  the  ECI  frame  to  coordinates  in  the  RTN  frame,  the  argument  of  perigee,  the 
right  ascension  of  the  ascending  node,  the  true  anomaly,  and  the  inclination  of  the  orbit 
must  be  known.  The  only  one  of  these  elements  that  varies  over  the  course  of  the  orbit  is 
the  true  anomaly  and  it  can  be  calculated  by 

u  =  (rem(T  I  P))2tt  (3.15) 


where: 

•  v  is  the  true  anomaly 
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•  T  is  the  time  passed  since  beginning  of  simulation 

•  P  is  the  orbital  period 

This  process  will  determine  the  true  anomaly  of  the  given  satellite  at  the  desired  instant  in 
time.  Therefore  the  position  and  velocity  vectors  in  the  ECI  frame  can  be  converted  to 
the  RTN  frame  by  the  following  transformations 


rRTN  ^ot^i^n^i,reci 


VRTN  ^co^i^n^uVeci 


(3.16) 

(3.17) 


where  the  rotation  matrices  are  as  follows: 


R  = 


cos(<z>)  sin(<z>)  0 

-sin(<x>)  cos(<x>)  0 

0  0  1 


(3.18) 


R:  = 


10  0 
0  cos(z')  sin(z') 

0  -sin(z')  cos(z') 


(3.19) 


cos(Q)  sin(Q)  0 
-sin(Q)  cos(Q)  0 
0  0  1 


(3.20) 
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R. 


(3.21) 


cos(y)  sin(f;)  0 

-sin(u)  cos  (u)  0 
0  0  1 


In-Plane  Formation 

For  the  in-plane  formation  the  deputy  satellite  has  all  of  the  same  initial 
conditions  as  the  control  satellite  except  that  it  is  a  pre-determined  distance  (1  km  for  this 
study)  behind  the  control  satellite  in  the  same  orbital  plane.  So  the  following  constraints, 
expressed  in  RTN  frame  coordinates,  must  be  met  for  the  satellites  to  meet  the  in-plane 
formation  criteria: 


x{t )  =  0 

(3.22) 

y(t)  =  v0 

(3.23) 

z(t )  =  0 

(3.24) 

The  argument  of  latitude  separation  can  be  detennined  by  inputting  the  desired  separation 
distance  between  the  master  and  deputy  satellites  into  the  following  equation 


ud  =  s/r 


(3.25) 


where: 

•  ud  is  the  argument  of  latitude  of  the  deputy  satellite 

•  s  is  the  separation  distance  between  the  deputy  and  control 
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•  r  is  the  orbital  radius 


In-Track  Formation 

For  the  in- track  experiments  the  deputy  satellite  not  only  has  a  different  argument 
of  latitude,  but  it  also  has  a  different  right  ascension  of  the  ascending  node  such  that  its 
ground  track  matches  that  of  the  control  satellite’s.  Hill’s  equations  for  this  formation  are 
very  similar  to  the  in-plane  formation  except  that  a  cross-track  oscillation  represents  the 
difference  in  right  ascension  of  the  ascending  node.  Hill’s  equations  for  this  type  of 
formation  can  be  expressed  in  the  RTN  frame  as: 


x(t)  =  0 

(3.26) 

y(0  =  To 

(3.27) 

z(t)  =  -(coe  /  n)y0  sin  /'  cos  nt 

(3.28) 

To  calculate  the  correct  right  ascension  of  the  ascending  node  to  keep  the  deputy  satellite 
following  the  control  satellite’s  ground  track  the  following  equation  was  used 

CLd={(2n-ud)l{27t))TcoE  (3.29) 


where: 

•  Cld  is  the  right  ascension  of  the  ascending  node  for  the  deputy 

•  ud  is  the  argument  of  latitude  of  the  deputy  satellite 

•  T  is  period  of  the  orbit 
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•  coE  is  the  rotational  rate  of  the  Earth 


Circular  Formation 

As  previously  discussed  the  circular  formation  is  designed  such  that  the  satellites 
will  maintain  a  constant  distance  from  each  other.  By  examining  Hill’s  equations  the 
relation  between  the  initial  conditions  must  meet  the  following  constraint  equation  where 
r  is  the  radius  of  the  fonnation  from  the  reference  point  at  the  control  satellite: 


2,2,2  2 
x  4-  y  +  z  -r 

(3.30) 

The  geometric  approach  takes  advantage  of  the  fact  that  the  relative  motion  in  the 

radial/tangential  planes  (x/y)  is  fixed  in  eccentricity.  By  substituting  the  constraints  into 

Hill’s  equations  the  following  relations  are  found. 

To  =  -2»*0 

(3.31) 

y0  =  2x0/n 

(3.32) 

o 

+i 

ii 

o 

N 

(3.33) 

z0  =  ±V3i0 

(3.34) 

The  first  two  equations  set  the  along  track  (tangential)  drift  and  offset  to  zero.  The 
second  two  equations  must  be  of  the  same  sign  for  any  particular  case.  These  constraints 
demonstrate  that  there  are  two  planes  in  which  a  circular  formation  is  possible. 
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Initial  Conditions  (Circular  Orbit) 


For  the  test  scenario  the  design  will  mimic  the  design  presented  by  Sabol,  et  al 
(2001 :  272)  in  their  paper  and  recounted  earlier  in  the  literature  review.  The  reference 
orbit  is  circular  and  inclined  at  90  degrees  with  the  satellite  at  the  equator  (ascending 
node).  The  deputy  satellite  is  also  at  the  equator,  but  at  its  ascending  node  and  in  a 
slightly  elliptical  orbit.  The  two  satellites  are  also  offset  by  a  small  amount  of  right 
ascension  of  the  ascending  node.  Given  this  type  of  formation  where  both  satellites  have 
the  same  inclination,  it  also  helps  to  diminish  the  disruptive  effects  of  the  Earth’s 
oblateness. 

In  order  to  ensure  that  both  satellites  in  the  circular  formation  have  the  same 
period  (and  thus  will  not  quickly  drift  apart)  it  is  essential  to  the  formation  stability  that 
all  satellites  in  the  formation  have  the  same  orbital  energy.  Therefore  to  set  the  initial 
conditions  of  the  control  satellite  and  the  deputy  satellite  use  the  energy  equation 

E  =  —v2  (3.35) 

2  r 

and  set  the  energy  of  the  control  satellite  equal  to  the  energy  of  the  deputy  satellite.  For 
the  initial  conditions  in  our  test  scenario,  the  difference  in  position  is  due  to  differences  in 
the  radial  and  normal  direction  and  the  difference  in  velocity  is  solely  due  to  a  delta  in  the 
tangential  direction.  Setting  the  energy  (and  consequently  the  period)  of  the  two  orbits  to 
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equal  values  ensures  that  the  satellites  will  remain  in  the  same  vicinity  and  that  tiny 
differences  in  the  orbital  energy  will  not  interfere  with  the  simulation. 


Controller  (Ogata,  1970:  156-157) 

An  automated  controller  compares  the  actual  value  of  the  system  to  the  desired 
value,  determines  the  difference,  and  produces  a  modifying  control  that  will  minimize  or 
eliminate  the  difference.  In  this  experiment  a  proportional  and  integral  controller  will  be 
used  help  keep  the  clustered  satellites  in  formation.  For  a  controller  with  a  proportional 
control  action  the  relationship  between  the  output  of  the  controller  m(t)  and  the  actuating 
error  signal  e(t)  is: 


m{t)  =  Kpe(t ) 


(3.36) 


where  K  is  tenned  the  proportional  sensitivity  or  the  gain.  This  type  of  system  can  be 

used  in  many  different  environments,  but  whatever  the  actual  mechanism  and  operating 
power,  the  proportional  controller  is  basically  an  amplifier  with  an  adjustable  gain. 

In  addition  to  the  proportional  control  action,  an  integral  control  action  is  also 
employed  such  that  the  control  of  the  satellite  will  become  more  stable.  In  a  controller 
with  integral  control  action,  the  rate  of  change  of  the  controller  output  m(t)  is 
proportional  to  the  actuating  error  signal  e(t): 


dm(t) 

dt 


K,e(t) 


(3.37) 
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For  the  simulations  done  in  this  experiment  the  following  values  were  chosen  to  obtain 
fairly  well  controlled  results: 


Kp  =150000 

(3.38) 

Kt  =  20000 

(3.39) 

This  controller  is  then  used  to  continuously  alter  the  cross  sectional  area  of  the  deputy 
satellite  so  that  it  remains  consistent  with  the  cross  sectional  area  of  the  control  satellite. 
This  controlling  action  should  mitigate  the  differential  drag  that  was  originally  present 
and  allow  the  satellites  to  remain  relatively  close. 

This  type  of  control  (proportional-integral)  is  typically  used  for  linear  systems.  In 
this  case,  it  has  been  applied  to  a  non-linear  system  and  thus  stability  cannot  be  assumed. 
During  the  scenarios  tested  for  this  thesis  the  controller  performed  well  despite  operating 
outside  of  its  intended  environment.  If  this  technique  is  pursued  in  the  future,  it  may  be 
necessary  and  prudent  to  perform  a  more  rigorous  analysis. 
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IV.  RESULTS 


Chapter  Overview 

This  section  is  broken  down  into  two  primary  parts.  The  first  part  deals  with  how 
drag  affects  satellites  in  each  of  the  three  discussed  formations  (in-plane,  in-track,  and 
circular).  Plots  and  analysis  will  be  included  that  illustrate  how  the  formation  is  reacting 
over  a  period  of  time.  Generally  the  analysis  and  plots  will  be  over  a  ten  hour  period; 
however  certain  configurations  will  require  longer  observation  periods  in  order  to  see  the 
disturbing  effects.  The  second  part  of  this  chapter  will  focus  on  using  a  controller  to 
alleviate  the  drift  that  the  differential  drag  is  causing  between  the  satellites  in  the  cluster. 

Un-Controlled  Formations 

In-Plane  Formation 

The  first  formation  examined  in  this  research  is  the  in-plane  formation.  As 
expected,  when  there  is  a  differential  drag  force  present  between  the  two  satellites  in  an 
in-plane  fonnation  they  begin  to  move  apart  and  will  continue  to  move  further  and  further 
apart  along  an  ever-increasing  curve.  The  following  graphs  show  the  separation  between 
the  two  satellites  in  the  radial,  tangential  and  normal  directions  as  well  as  the  total 
separation. 
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x  103  Radial  Direction  In-Plane 


Figure  4.  In-plane  radial  separation 


The  radial  direction  separation  in  the  in-plane  formation  begins  to  slowly  grow  as  the 
differential  drag  between  the  control  and  deputy  satellites  causes  the  control  satellite  to 
speed  up  and  enter  a  slightly  lower  orbit.  Figure  4  shows  that  after  ten  hours  the  radial 
separation  is  6.5  meters  and  growing  at  a  steady  rate.  This  separation  is  in  the  radial 
direction  and  means  that  the  control  satellite  now  has  a  slightly  smaller  semi-major  axis 
and  thus  a  different  period.  The  separation  in  this  direction  will  continue  to  grow  linearly 
unless  some  type  of  contrary  force  can  be  applied  or  the  differential  drag  can  be 
equalized. 
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Tangential  Direction  In-Plane 


Figure  5.  In-plane  tangential  separation 

Originally  the  two  satellites  are  one  kilometer  apart  in  the  tangential  direction  in  the  in¬ 
plane  formation.  However  because  the  satellite  in  the  formation  lead  (the  control 
satellite)  has  a  slightly  larger  cross  sectional  area  due  to  an  anomaly  on  the  deputy 
satellite  it  will  experience  a  larger  (negative)  acceleration  due  to  drag.  This  will  cause  the 
control  satellite’s  semi-major  axis  to  decrease  and  force  it  to  speed  up.  As  the  control 
satellite  continues  to  accelerate  away  from  the  deputy  the  tangential  separation  will 
continue  growing  at  an  ever-increasing  rate. 
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Normal  Direction  In-Plane 
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Figure  6.  In-plane  normal  separation 

The  in-plane  formation  calls  for  two  satellites  to  be  in  the  same  plane  and  separated  by 
mean  anomaly.  Their  separation  in  the  normal  direction  should  remain  at  zero  since  drag 
acts  solely  in  the  orbital  plane  and  no  other  perturbing  forces  were  considered.  Figure  6 
demonstrates  that  the  satellites  do  maintain  their  separation  of  zero  in  the  normal 
direction. 
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Separation  Between  Satellites(ln-Plane) 


Figure  7.  Total  separation  between  satellites  (in-plane) 


Because  the  two  satellites  occupy  the  same  orbital  plane,  their  separation  distance  will  be 
solely  in  the  radial  and  tangential  directions  (as  shown  in  Figure  4,  Figure  5,  and,  Figure 
6).  In  addition,  as  drag  begins  to  affect  the  control  satellite  more  than  the  deputy  satellite, 
the  control  satellite  will  begin  to  speed  up,  decrease  orbital  altitude,  and  move  away  from 
the  deputy  satellite  at  a  faster  and  faster  rate. 
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In-Track  Formation 


The  in-track  formation  radial  direction  separation  begins  to  slowly  grow  as  the 
drag  on  the  control  satellite  slowly  pulls  it  closer  to  Earth.  Figure  8  shows  that  after  ten 
hours  the  radial  separation  is  6.5  meters.  This  separation  is  in  the  radial  direction  and 
means  that  the  control  satellite  now  has  a  slightly  smaller  semi -major  axis  and  thus  a 
different  period.  The  separation  in  this  direction  will  continue  to  grow  linearly  unless 
some  type  of  contrary  force  can  be  applied  or  the  differential  drag  can  be  equalized. 


x103  Radial  Direction  In-Track 


Figure  8.  In-track  radial  separation 
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Tangential  Direction  In-Track 


Figure  9.  In- track  tangential  separation 

While  the  radial  separation  is  slowly  growing,  the  tangential  separation  is  also  growing  at 
a  similar  pace  for  the  first  ten  hours.  However  because  the  tangential  separation  is 
growing  along  a  quadratic  curve  it  will  soon  become  the  dominant  direction  of 
separation.  As  the  periods  of  the  two  satellites  diverge,  the  satellite  with  a  larger  cross 
sectional  area  (control  satellite  in  this  case)  will  begin  to  speed  up  and  pull  farther  and 
farther  away  from  the  deputy  satellite.  The  control  satellite  will  continue  to  accelerate 
away  from  the  deputy  satellite  and  the  tangential  separation  will  continue  to  increase  at 
an  ever-increasing  rate. 
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Normal  Direction  In-Track 


Figure  10.  In- track  normal  separation 


Because  the  two  satellites  have  slightly  different  ascending  nodes,  their  separation  in  the 
normal  direction  is  periodic,  with  a  period  equal  to  that  of  the  orbits.  This  is  illustrated  in 
Figure  10.  The  magnitude  of  the  separation  will  vary  depending  on  the  altitude  of  the 
orbits.  Fligher  satellite  orbits  will  increase  the  magnitude  of  the  separation  (assuming 
their  initial  separations  are  equal  in  each  case). 
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In-Track  3-D 


Figure  11.  In-track  3-D  (100  hours) 


Figure  1 1  demonstrates  how  the  in-track  formation  changes  over  100  hours.  As 
expected,  the  tangential  separation  becomes  the  dominant  direction  of  separation. 
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Circular  Formation 


For  the  circular  formation  the  orbits  were  also  initially  propagated  for  10  hours. 
However,  some  of  the  motion  was  not  easily  discernable  over  the  ten  hour  period  so 
additional  plots  are  included  with  different  parameters  to  show  the  motion  of  the  satellites 
in  the  circular  formation. 


Radial  Direction  Circular 


Figure  12.  Circular  radial  separation 


As  previously  described,  in  the  circular  formation  both  the  deputy  and  control  satellite 
have  the  same  period  but  the  deputy  satellite  has  a  slightly  eccentric  orbit  that  causes  it  to 
oscillate  around  the  control  satellite  in  regards  to  the  radial  direction.  The  chosen  radial 
separation  at  the  beginning  of  the  simulation  was  0.5  km.  Because  the  deputy  satellite 
was  at  apogee,  that  separation  distance  should  also  be  the  maximum  separation. 
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Given  the  scale  of  Figure  12,  and  the  small  difference  in  cross  sectional  areas 
between  the  control  and  deputy,  it  is  difficult  to  discern  any  change  in  maximum 
separation  between  the  two  satellites.  To  better  illustrate  how  the  formation  will  change 
over  time  the  cross  sectional  area  difference  is  exaggerated  in  the  case  shown  in  Figure 
13  to  illustrate  how  the  formation  will  change  over  time  as  the  differential  drag  affects 
the  formation. 


Radial  Direction  Circular 


Figure  13.  Circular  radial  separation  (exaggerated) 


The  control  satellite’s  altitude  drops  faster  than  the  deputy’s  altitude  because  of  the 
greater  deceleration  due  to  drag  that  it  experiences.  As  the  semi-major  axes  change  the 
radial  separation  between  the  two  satellites  will  increase. 
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Radial  Direction  Circular 
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Figure  14.  Circular  radial  separation  (100  hours) 


It  is  not  obvious  from  Figure  14  that  the  radial  separation  is  changing  over  the  100  hour 
time  period  because  of  the  large  scale  that  must  be  used  to  show  the  radial  separation 
distance  oscillating.  However  in  Figure  15  the  top  edge  of  the  plot  is  amplified  so  that  it 
becomes  apparent  that  the  separation  distance  in  the  radial  direction  is  increasing. 
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Radial  Direction  Circular 


0.51 
0.508 
0.506 
0.504 

c„  0.502 

U— 

CD 

|  0.5 

_o 

2  0.498 
0.496 
0.494 
0.492 
0.49 

10  20  30  40  50  60  70  80  90  100 

Time  (hours) 

Figure  15.  Circular  radial  separation  zoomed  (100  hours) 

The  radial  separation  is  not  growing  quickly,  but  it  is  growing.  Over  the  100  hour 
propagation  time  the  radial  separation  grew  by  slightly  less  than  two  meters.  This  change 
is  expected  because  of  the  differential  drag  forces  that  the  two  satellites  are  experiencing. 


51 


Tangential  Direction  Circular 


Figure  16.  Circular  tangential  separation 


At  the  start  of  the  simulation  both  the  control  and  deputy  satellite  have  the  same  semi¬ 
major  axis.  The  control  satellite  is  in  a  circular  orbit  and  the  deputy  is  in  slightly 
elliptical  orbit  so  their  tangential  separation  will  oscillate  depending  on  where  the  two 
satellites  are  located  in  their  respected  orbits.  Their  maximum  separation  will  occur 
when  the  deputy  satellite  is  either  at  apogee  or  perigee. 

Given  the  scale  of  Figure  16,  and  the  small  difference  in  cross  sectional  areas 
between  the  control  and  deputy,  it  is  difficult  to  discern  any  change  in  maximum 
tangential  separation  between  the  two  satellites.  To  better  illustrate  the  changes  to  the 
formation  over  time  a  graph  with  exaggerated  parameters  is  shown  in  Figure  17. 
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Tangential  Direction  Circular 


Figure  17.  Circular  tangential  separation  (exaggerated) 


The  tangential  separation  in  the  circular  formation  begins  to  grow  in  a  similar  manner  to 
the  tangential  separation  in  the  in-plane  and  in-track  formations.  However  because  of  the 
elliptical  nature  of  the  deputy’s  orbit  there  is  an  oscillation  in  the  separation  distance. 

The  control  satellite  continues  to  drop  in  altitude  due  to  the  drag  deceleration  and  as  it 
drops  it  accelerates  away  from  the  deputy  satellite.  Therefore  the  faster  the  control 
satellite  goes,  the  farther  back  the  deputy  will  fall  and  the  tangential  separation  will 
continue  to  grow  along  a  quadratic  curve  unless  the  differential  drag  can  be  corrected. 
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Figure  18.  Circular  tangential  separation  (100  hours) 

Even  with  the  initial  cross  sectional  areas  (control  =  2  meters  squared  and  deputy  =1.5 
meters  squared)  a  major  change  in  the  tangential  separation  is  evident  if  the  orbits  are 
propagated  out  over  100  hours.  The  distance  continues  to  oscillate  around  a  mean  value, 
but  as  more  time  passes  the  mean  value  of  the  oscillation  slowly  drifts  away  from  zero. 
After  100  hours  the  center  of  the  oscillation  has  drifted  nearly  500  meters  away  from  the 
starting  position.  The  tangential  separation  quickly  becomes  the  dominant  direction  of 
separation  in  the  circular  formation. 
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Normal  Direction  Circular 


Figure  19.  Circular  normal  separation 


Because  the  two  satellites  have  slightly  different  ascending  nodes,  their  separation  in  the 
normal  direction  is  periodic  with  a  period  equal  to  the  orbital  period.  This  is  illustrated  in 
Figure  19.  Because  the  drag  force  acts  in  a  direction  opposite  the  velocity  vector,  and 
both  the  control  and  deputy  satellites  are  in  a  polar  orbit  there  should  be  no  appreciable 
change  to  the  normal  separation  over  the  life  of  the  formation. 
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Kilometers 


Normal  Direction  Circular 


Figure  20.  Circular  tangential  separation  (exaggerated) 

Looking  at  the  separation  in  the  normal  direction,  Figure  20  confirms  the  belief  that 
differential  drag  will  not  appreciably  affect  the  separation  in  the  normal  direction. 
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Circular  3-D 


Figure  21.  Circular  tangential  separation  with  small  drag  differences  (10  hours) 

Figure  2 1  is  a  three  dimensional  plot  of  the  motion  of  the  two  satellites  in  a  circular 
formation.  The  control  satellite  has  a  cross  sectional  area  of  2  meters  squared  and  the 
deputy  has  a  cross  sectional  area  of  1.5  meters  squared.  The  circular  formation  appears 
to  maintain  its  integrity  pretty  well  for  the  first  ten  hours  of  the  simulation.  Figure  22  and 
Figure  23  show  the  motion  of  the  formation  over  100  hours  using  the  same  cross 
sectional  areas  for  the  control  and  deputy  satellites. 
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Circular  3-D 


Figure  22.  Circular  3-D  separation  with  small  drag  differences  (100  hours) 


Figure  22  is  a  three  dimensional  plot  of  the  motion  between  the  two  satellites  in  the 
circular  formation  over  a  100  hour  time  period.  Despite  the  differential  drag  situation  the 
formation  still  maintains  a  somewhat  circular  formation;  however  it  is  obvious  that  the 
formation  is  beginning  to  degrade. 
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Circular  3-D 


Figure  23.  Circular  3-D  separation  with  small  drag  differences  rotated  (100  hours) 

Figure  23  shows  the  same  case  as  Figure  22  from  a  different  angle.  In  this  figure  (Figure 
23)  it  is  clear  that  the  formation  degrades  primarily  in  the  tangential  direction.  Over  time, 
the  tangential  separation  will  continue  to  dominate  the  total  separation  distance. 
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Radial  Direction  (Kilometers) 


Circular  3-D 


Figure  24.  Circular  3-D  separation  (exaggerated,  10  hours) 

Earlier  in  this  section  a  few  plots  were  displayed  where  the  control  and  deputy  satellites 
had  exaggerated  differences  in  their  cross  sectional  areas  (500  meters  squared  vs.  1.5 
meters  squared).  Figure  24  shows  a  three  dimensional  combination  of  those  graphs. 
Once  again  it  is  clear  that,  as  the  control  satellite  loses  altitude  and  speeds  up,  it  will 
cause  ever  increasing  separations  in  the  tangential  direction. 
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Controlled  Formations 


The  previous  section  displayed  graphs  and  analysis  describing  how  different 
formations  reacted  to  satellites  in  the  formation  experiencing  different  drag  forces.  In  all 
cases  it  was  obvious  that  over  a  period  of  time  the  satellites  fell  out  of  formation  as  their 
positions  continued  to  diverge.  This  presents  a  problem  for  fonnations  whose 
deployment  does  not  go  perfectly. 

If  one  satellite  in  the  fonnation  has  a  solar  array  that  does  not  completely  deploy, 
or  a  piece  of  thermal  blanketing  comes  loose,  or  an  antenna  is  bent  slightly,  or  any 
number  of  other  things  happen  that  change  the  physical  configuration  of  a  satellite  then 
one  (or  more  than  one)  of  the  satellites  in  the  fonnation  will  have  a  slightly  different 
cross-sectional  area  than  the  normal  and  it  will  experience  a  different  acceleration  due  to 
drag  This  differential  drag  will  eventually  disrupt  the  configuration  of  the  cluster  and 
possibly  jeopardize  the  mission  of  the  clustered  satellites. 

However,  if  this  slight  difference  in  cross  sectional  area  can  be  corrected  then  the 
two  satellites  will  not  experience  the  differential  drag  forces  and  they  will  remain  in  the 
formation  (subject  to  other  perturbations).  In  order  to  accomplish  this  effect  a 
proportional  and  integral  controller  was  introduced  into  the  code.  The  controller 
compares  the  specific  mechanical  energy  of  the  control’s  and  deputy’s  orbit  and  then 
assigns  a  cross  sectional  area  to  the  deputy  satellite  that  is  proportional  to  the  energy 
difference.  This  type  of  controller  could  be  implemented  on  a  satellite  by  adding  drag 
plates  that  can  be  shifted  to  different  angles  in  relation  to  the  velocity  vector  to  change 
the  satellite’s  cross  sectional  area.  Unless  otherwise  specified  all  of  the  following 
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graphs  have  been  generated  using  the  following  controller  gains  for  the  proportional  and 


integral  controllers: 


Kp  =150000 

(4.1) 

K.  = 20000 

(4.2) 

In  addition,  a  controller  with  much  smaller  gains  will  be  characterized  on  the  in-plane 
scenarios  to  determine  if  a  less  capable  controller  can  still  maintain  the  desired 
configurations.  This  controller  has  proportional  and  integral  gains  of: 


Kp  =  15 

(4.3) 

K,=  2 

(4.4) 

When  the  controller  gains  specified  by  equations  (4.3)  and  (4.4)  is  used  to  generate 
graphs  the  proportional  and  integral  control  values  will  be  specified  on  the  graph. 


Controlled  In-Plane  Formation 

The  controller  discussed  above  and  in  the  methodology  chapter  was  used  to  change  the 
cross  sectional  area  of  the  deputy  satellite  in  the  following  graphs.  The  initial  values  for 
cross  sectional  areas  ( 1 .5  meters  squared  for  the  deputy  and  2  meters  squared  for  the 
control)  are  displayed  on  each  graph. 
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x  103  Radial  Direction  In-Plane 


Figure  25.  In-plane  radial  separation  with  controller  (100  hours) 


In  the  previous  uncontrolled  example  for  the  in-plane  formation  there  was  an  obvious 
divergence  between  the  control  and  deputy  satellite  after  10  hours.  In  this  case  the 
scenario  was  run  for  100  hours  and  the  controller  demonstrates  the  ability  to  correct  for 
the  drag  differences  and  keep  the  formation  close  together  in  the  radial  direction. 
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Tangential  Direction  In-Plane 


Figure  26.  In-plane  tangential  separation  with  controller  (100  hours) 


For  the  uncontrolled  in-plane  formation  the  tangential  separation  began  to  quickly  grow 
and  continued  growing  along  a  quadratic  curve.  However,  with  the  addition  of  the 
controller,  the  tangential  separation  stays  within  centimeters  of  the  starting  separation  (1 
km). 
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Normal  Direction  In-Plane 
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Figure  27.  In-plane  radial  separation  with  controller  (100  hours) 

The  separation  in  the  normal  direction  did  not  increase  in  the  uncontrolled  scenario  so 
there  is  no  reason  to  believe  that  it  will  change  in  the  controlled  scenario.  The  above 
graph  demonstrates  that  the  normal  separation  remains  right  at  zero  as  expected. 
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Tangential  vs.  Radial  Separation  In-Plane 
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Figure  28.  In-plane  tangential  vs.  radial  separation  with  controller  (100  hours) 

As  expected  when  the  tangential  and  radial  separation  plots  are  combined,  it  is  easy  to  see 
that  the  formation  is  stable  when  the  controller  is  functioning  properly.  A  three- 
dimensional  plot  is  not  necessary  because  the  separation  in  the  normal  direction  is  zero 
for  all  intents  and  purposes. 
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Cross  Sectional  Area  vs.  Time 


Figure  29.  Cross  sectional  area  vs.  time  maximum  control  (100  hours) 

Figure  29  shows  how  the  controller  alters  the  deputy  satellite’s  cross  sectional  area  such 
that  the  orbits  of  the  control  and  deputy  satellites  have  the  same  energy.  Because  this 
particular  controller  is  capable  of  producing  large  gains,  the  deputy  satellite  changes 
cross  sectional  area  at  a  rapid  pace  and  the  separation  distances  can  be  kept  to  a 
minimum.  In  addition,  as  this  controller  continues  to  operate  it  will  require  less  drastic 
cross  sectional  area  changes  and  continue  to  oscillate  closer  to  the  cross  sectional  area  of 
the  control  satellite.  The  controller  used  for  this  study  only  used  proportional  and 
integral  control;  however,  the  addition  of  a  derivative  control  and  further  tuning  of  the 
controller  gains  would  both  help  dampen  the  oscillation  and  allow  the  cross  sectional  area 
to  settle  down  much  more  quickly  (Ogata,  1970:  157).  The  derivative  control  action  has 
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an  anticipatory  character  and  is  only  effective  during  transition  periods  (Ogata,  1970: 
158). 

The  next  step  for  examining  the  controlled  in-plane  formation  involves  examining 
how  the  formation  will  react  to  a  controller  that  is  not  capable  of  generating  large  enough 
gains  to  immediately  correct  for  the  differential  drag.  The  following  graphs  were 
propagated  over  longer  periods  of  time  so  that  the  affect  of  the  controller  could  be  better 
characterized.  The  cross  sectional  areas  displayed  on  the  graphs  are  the  initial  cross 
sectional  areas,  however  because  the  controller  adjust  the  cross  sectional  areas  it  will 
change  constantly  throughout  the  simulation.  The  change  in  cross  sectional  area  is 
presented  in  a  later  graph. 


x  10-4  Radial  Direction  In-Plane 


Figure  30.  In-plane  radial  separation  with  minimum  control  (500  hours) 
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As  the  control  satellite’s  orbit  lowers  the  controller  must  adjust  the  cross  sectional  area  of 
the  deputy  satellite  such  that  both  orbits  have  the  same  energy.  To  do  this  the  controller 
increases  the  cross  sectional  area  of  the  deputy  satellite.  The  less  powerful  controller 
employed  in  this  case  cannot  adjust  the  cross  sectional  area  immediately  (as  the  previous 
controller  did).  As  the  controller  adjusts  the  cross  sectional  area  it  puts  the  deputy 
satellite  into  a  slightly  elliptical  orbit  so  that  the  two  satellites  will  maintain  the  same 
orbital  energy.  This  causes  the  radial  separation  of  the  two  satellites  to  vary  periodically 
over  a  period  of  -200  hours. 


Tangential  Direction  In-Plane 


Figure  3 1 .  In-plane  tangential  separation  with  minimum  control  (500  hours) 
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The  tangential  separation  between  the  two  satellites  follows  a  plot  very  similar  to  the 
radial  separation  plot.  It  is  slowly  increasing  as  the  controller  increases  the  cross 
sectional  area  and  then  slowly  moves  back  toward  the  nominal  separation. 
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Figure  32.  In-plane  normal  separation  with  minimum  control  (500  hours) 

The  normal  separation  continues  to  be  a  non-issue  and  for  all  intents  and  purposes  it  can 
be  treated  as  zero. 
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3.5 


Cross  Sectional  Area  vs.  Time 


Figure  33.  Cross  sectional  area  vs.  time  with  minimum  control  (500  hours) 


At  the  initiation  of  the  scenario,  the  control  satellite  has  a  larger  cross  sectional  area  than 
the  deputy  satellite  and  thus  the  controller  begins  to  increase  the  cross  sectional  area  of 
the  deputy  satellite  in  an  attempt  to  equalize  the  orbital  energies.  By  the  time  the 
controller  has  equalized  the  cross  sectional  areas  of  the  two  satellites  the  control  satellite 
is  already  in  a  lower  orbit.  Hence,  the  controller  continues  increasing  the  cross  sectional 
area  of  the  deputy  satellite  as  it  works  toward  equalizing  the  orbital  energies.  When  the 
deputy  satellite’s  cross  sectional  area  reaches  approximately  2.5  meters  squared  the 
orbital  energy  of  the  two  orbits  are  equal;  however,  the  controller  cannot  reduce  the  cross 
sectional  area  of  the  deputy  fast  enough  to  keep  the  energies  equal.  Because  the  deputy 
satellite  has  a  larger  cross  sectional  area  than  the  control  satellite  it  will  fall  into  a  slightly 
lower  orbit  and  the  controller  will  begin  to  reduce  the  cross  sectional  area  of  the  deputy 
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until  it  reaches  1 .5  meters  squared  again.  During  this  time  the  tangential  separation  is 
being  reduced  and  the  two  satellites  are  moving  back  to  the  nominal  separation  of  1  km. 
However  when  the  satellites  reach  the  nominal  separation  of  1  km  the  cross  sectional 
areas  are  also  back  at  their  initial  values  and  the  entire  process  starts  over.  Thus  the 
periodic  motion  shown  in  the  radial  and  tangential  separation  plots  is  also  mirrored  in  the 
cross  sectional  area  versus  time  plot  that  is  displayed  above. 

This  periodic  motion  could  possibly  be  mitigated  with  the  addition  of  a  derivative 
control  action  (sometimes  called  rate  control),  or  by  better  tuning  the  K  and  KI  gains. 

The  derivative  control  action  has  an  anticipatory  character  and  is  effective  during 
transient  periods;  however  the  derivative  control  can  never  be  used  alone  and  must  be 
implemented  in  conjunction  with  a  proportional  or  a  proportional  plus  integral  control 
action  (Ogata,  1970:  157). 

Controlled  In-Track  Formation 

The  in-track  formation  demonstrated  a  diverging  tendency  similar  to  the  in-plane 
formation.  When  the  controller  was  added  to  the  scenario  it  corrected  the  divergence  but 
the  results  were  slightly  different.  The  in-track  scenarios  were  all  propagated  using  the 
original  controller  with  large  gains. 
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Radial  Direction  In-Track 


x  10'3 


Figure  34.  In-track  radial  separation  with  controller  (100  hours) 


In  the  uncontrolled  scenario  the  radial  separation  changed  by  about  6.5  meters  over  a  10 
hour  time  period  shown.  When  the  controller  was  employed,  it  reduced  the  maximum 
separation  to  less  than  one  centimeter  as  shown  by  the  above  graph.  If  the  controller  is 
allowed  to  continue  working,  the  maximum  separation  will  continue  to  decrease  over 
time. 
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Radial  Direction  In-Track 


x  10" 


Figure  35.  In- track  radial  separation  with  controller  (100  hours  zoom  in) 


When  the  in-track  scenario  was  propagated  over  100  hours,  it  becomes  obvious  that  the 
radial  separation  will  continue  to  decrease  as  the  controller  is  allowed  to  operate  and 
slowly  reduce  the  oscillation  to  a  near  zero  value. 
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Tangential  Direction  In-Track 


Figure  36.  In- track  tangential  separation  with  controller  (100  hours) 


In  the  uncontrolled  scenario  the  tangential  separation  began  to  grow  at  an  ever  increasing 
pace  after  a  few  hours;  however,  when  the  controller  was  implemented,  the  separation 
could  be  contained  to  a  negligible  value. 
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Tangential  Direction  In-Track 


Figure  37.  In-track  tangential  separation  with  controller  (100  hours  zoom  in) 

Propagating  the  in-track  scenario  over  100  hours  and  zooming  in  to  the  area  of  interest 
yields  Figure  37.  This  graph  is  identical  to  Figure  36  only  with  a  smaller  range  of  values 
along  the  Y-axis.  Clearly,  as  the  controller  continues  to  operate,  the  tangential  separation 
will  settle  down  closer  and  closer  to  a  one  kilometer  separation. 
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Normal  Direction  In-Track 


Figure  38.  In-track  normal  separation  with  controller  (10  hours) 

For  the  in-track  formation  the  only  differences  between  the  orbits  of  the  deputy  and 
control  satellites  are  in  the  mean  anomaly  and  the  right  ascension  of  the  ascending  node. 
The  difference  in  right  ascension  of  ascending  node  will  cause  the  normal  separation  to 
oscillate  around  zero  with  the  largest  differences  being  when  the  satellites  are  near  the 
equator. 

A  figure  similar  to  Figure  28  could  be  displayed  for  the  in-track  scenario; 
however,  the  figure  looks  identical  to  that  figure  and  thus  has  not  been  included. 
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Controlled  Circular  Formation 


The  circular  formation  was  tested  with  exaggerated  differences  between  the  cross- 
sectional  areas  of  the  two  satellites  to  show  how  the  formation  would  degrade 
propagating  the  scenario  over  a  reasonable  amount  of  time.  While  it  is  not  reasonable  to 
imagine  two  satellites  identical  upon  launch  having  such  a  wide  disparity  upon 
deployment,  it  does  illustrate  the  trends  of  the  formation  quite  well  in  an  easy  to  read  and 
understand  plot.  For  the  controlled  scenarios  the  huge  cross-sectional  area  differences 
will  continue  to  be  used  to  test  the  controller. 


Radial  Direction  Circular 


Figure  39.  Circular  radial  separation  with  controller  (10  hours) 


78 


For  the  uncontrolled  scenario  the  separation  began  to  expand  in  a  noticeable  manner  after 
the  plotted  ten  hour  period.  However  with  the  controller  employed  the  radial  separation 
maintains  the  constant  oscillation  that  is  expected  of  two  satellites  with  the  same  semi¬ 
major  axis  yet  slightly  different  eccentricities. 


Tangential  Direction  Circular 


Figure  40.  Circular  tangential  separation  with  controller  (10  hours) 


The  tangential  direction  of  separation  quickly  diverged  from  zero  when  tested  in  the  un¬ 
controlled  scenario,  but  as  expected  once  the  controller  was  added  to  the  scenario  it 
settled  down  and  oscillated  around  zero  without  diverging.  The  points  of  maximum 
separation  occur  near  the  north  and  south  poles  (which  is  consistent  with  the  design  of  the 
circular  formation  used  for  these  scenarios). 
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Normal  Direction  Circular 


Figure  41 .  Circular  normal  separation  with  controller  (10  hours) 


Once  again  no  major  changes  occurred  in  the  normal  direction  of  the  un-controlled 
scenario  so  no  changes  are  expected  in  the  controlled  scenarios.  Figure  41  shows  that  the 
separation  in  the  normal  direction  is  no  different  than  the  un-controlled  scenario  and  it 
continues  to  oscillate  around  zero  depending  upon  which  part  of  the  orbit  the  satellites 
occupy  at  the  current  time. 
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Figure  42.  Circular  3-D  separation  with  controller  (10  hours) 


Figure  42  is  a  three  dimensional  plot  of  the  controlled  circular  formation  with  the 
exaggerated  cross-sectional  area  differences.  In  the  un-controlled  scenario  the  tangential 
direction  of  separation  quickly  expanded  and  the  circular  formation  degraded.  The 
addition  of  the  controller  has  solved  that  problem  and  the  circular  formation  remains 
stable  while  the  controller  is  applied. 
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V.  Conclusions  and  Recommendations 


Conclusions 

The  objective  of  this  research  was  to  study  the  effects  of  drag  on  anomalous 
satellites  within  a  cluster  and  determine  how  quickly  the  cluster  formation  will 
deteriorate.  Obviously  this  question  depends  a  lot  on  the  drag  differential  present 
between  the  satellites.  This  study  examined  two  satellites  with  a  25%  difference  in  cross 
sectional  area  and  found  (as  expected)  that  the  fonnations  quickly  deteriorated  and  the 
primary  direction  of  separation  was  in  the  tangential  direction.  It  quickly  became 
obvious  that  none  of  the  studied  formations  would  be  able  to  maintain  their  fonnation 
without  some  type  of  controlling  mechanism. 

Of  the  three  formations  (in-plane,  in-track,  and  circular)  examined  in  this 
research,  the  satellites  in  the  circular  formation  maintained  proximity  for  the  longest 
period  of  time;  however,  this  research  did  not  include  the  effects  of  the  J2  perturbation 
which  would  greatly  affect  the  stability  of  the  circular  formation. 

Many  future  formations  will  require  precise  satellite  positioning  and  any  small 
differences  in  cross  sectional  areas  will  cause  the  satellites  to  diverge  from  the  designed 
formation.  Therefore,  if  some  type  of  controlling  mechanism  can  be  implemented  to 
alleviate  the  differential  drag  between  the  satellites,  the  formation  may  be  able  to  retain 
its  integrity  and  functionality  for  extended  periods  of  time. 

This  study  implemented  a  controller  that  adjusted  the  cross  sectional  areas  of  the 
satellites  such  that  the  energies  of  the  orbits  would  remain  equal.  If  this  type  of  controller 
can  be  implemented  with  large  enough  gains  then  the  differential  drag  can  be  corrected 
immediately  and  the  formation  will  remain  true  to  its  design.  However  if  a  controller 
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with  smaller  gains  is  used  for  the  controlling  actions  the  formations  will  degrade 
(depending  of  the  cross  sectional  area  difference)  before  the  controller  can  fully  correct 
the  problem.  Then  the  satellites  will  remain  in  close  proximity  but  the  integrity  of  the 
formation  design  will  not  be  maintained  and  the  total  separation  between  the  satellites 
will  vary  from  the  initial  separation  distance  to  an  arbitrary  maximum  distance  based  on 
the  drag  differential  and  the  controller  gain.  For  this  research,  the  arbitrary  distance  was 
20%  larger  than  the  initial  separation  and  the  separation  between  the  two  satellites  varied 
along  a  sine  curve  between  the  initial  separation  and  maximum  separation. 

At  some  point  satellite  clusters  will  become  imperative  to  the  advancement  of 
military  ISR  systems  and  at  that  point  a  controlling  system  that  can  compensate  for  drag 
differences  and  other  perturbations  will  also  become  necessary  to  ensure  the  success  of 
those  systems. 

Recommendations  for  Future  Research 

This  research  characterized  the  manner  in  which  the  separation  between  clustered 
satellites  in  different  fonnations  would  diverge  if  one  of  the  satellites  experienced  an 
anomaly  which  effected  its  cross  sectional  area.  In  addition,  a  proportional/integral 
controller  was  employed  to  compensate  for  the  different  cross-sectional  areas.  This  type 
of  controller  is  fairly  straight-forward  and  the  possibility  of  expanding  the  scope  and 
complexity  of  the  controller  would  be  interesting.  Adding  a  differential  control 
component  to  the  controller  and  fine  tuning  the  proportional  and  integral  gains  of  the 
controller  would  be  a  logical  step  forward.  Further  verification  of  the  controller  in  a  non¬ 
linear  environment  may  also  be  necessary.  Studying  the  dynamics  and  consumable 
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budgets  of  a  formation  that  is  forced  to  depend  on  this  style  of  controller  to  maintain  the 
functional  formation  would  also  be  a  logical  step  forward. 

Increasing  the  cross  sectional  area  to  mitigate  the  differential  drag  is  a  viable 
option  for  controlling  the  satellites  in  formation,  but  it  may  also  reduce  the  lifetime  of  the 
formation  (more  drag  =  shorter  lifetime).  Additional  research  could  be  done  that  studies 
the  impact  on  formation  lifetime.  Along  those  same  lines,  a  study  that  examined  how 
attitude  control  could  be  used  to  correct  for  differential  drag  would  be  very  interesting. 

Another  possible  research  area  is  the  adaptation  of  a  more  accurate  atmospheric 
model  into  the  code.  The  current  atmosphere  model  is  satisfactory  for  the  context; 
however,  results  obtained  from  a  higher  fidelity  model  would  be  important  in  the  long 
run.  Along  similar  lines,  a  study  that  implemented  the  full  scope  of  orbital  perturbations 
instead  of  just  two  body  motion  and  atmospheric  drag  would  provide  higher  fidelity 
results  and  may  uncover  some  different  issues  related  to  the  stability  and  behavior  of 
these  formation  designs. 

In  addition  to  the  three  types  of  formations  discussed  in  this  research,  there  are 
other  possible  formation  designs  that  could  be  explored.  Other  formation  types  will 
likely  require  extensive  formation  keeping.  Unless  those  hurdles  can  be  overcome  in  an 
efficient  manner,  there  may  not  be  a  need  to  develop  a  differential  drag  controller  to  be 
implemented  on  the  formation  satellites. 
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Appendix  A  -  Characterization  Code 


The  following  code  was  written  and  executed  using  MATLAB  version  7.0. 
function  []=orbit_propagator() 

%  Simulation  of  a  satellite  in  orbit  with  two-body  motion  only 
%  Initial  and  Final  Time  (seconds) 
tO  =  0; 

tf  =  36000;  %  in  seconds,  100  hrs  =  360000 
t_vec=t0:10:tf; 

%  Initial  Satellite  Orbital  Parameters 
%  Orbit  Altitude  (Kilometers) 

Sat_Alt  =  250; 

EarthRadius  =  6378.1; 

Sat_SemiMajor  =  Earth_Radius  +  Sat_Alt;  %  Compute  Semimajor  Axis 
%  Constants 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 

mu  =  3.98601e5; 

%  Gravitational  constant 
%  Expressed  in  N  mA2/kgA2 

G  =  6.67e-l  1; 

%  Mass  of  the  Earth  expressed  in  kg 
MassEarth  =  5.9742E24; 

%  Rotational  rate  of  Earth  (rad/s) 

Earth_rate  =  7.2722e-5;  %  corresponds  to  15  degrees/hour 
%  Orbit  Eccentricity 
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Sat_Ecc  =  0;  %  Eccentricity  assumed  zero 
%  Orbit  Inclination  (radians) 

Sat_Inc  =  pi/2; 

%  Right  Acsension  of  Ascending  Node  (radians) 

%  Angle  between  vernal  equinox  direction  and  point  at  which  the  satellite  passes 
%  the  equator  going  north 

Sat_RAAN  =  0; 

%  Argument  of  Perigee  (radians) 

%  Argument  of  Perigee  undefined  for  circular  orbit,  assume  zero  as  placeholder 
Sat_ArgPeri  =  0; 

%  Argument  of  Latitude  at  epoch  (radians) 

%  Angle  between  the  ascending  node  direction  and  the  satellite  position  vector 
%  at  the  start  time  of  the  sim 

Sat_ArgLat  =  0; 

%  Argument  of  Latitude  at  epoch  (radians) 

%  For  an  in-plane  satellite  that  is  approximately  1  km  behind  the 
%  control  satellite.  Calculation  2pi*radius  of  sat  orbit  = 

%  circumference  of  orbit 

%  circumference/360  =  1 15.68  km/degree 

%  1/1 15.68  =  .008644  degrees  for  a  one  km  separation. 

%  .008644*pi/180  =  .000150873 

PSat  ArgLat  =  2*pi  -  0.000150873; 

%  Angle  between  prime  meridian  and  vernal  equinox  direction  at  start  time  of  sim. 

%  This  is  an  arbitrary  assignment  unless  we  want  sim  time  to  have  some  real  meaning 
%  with  respect  to  a  true  astronomical  time  system 

EarthPos  =  0; 

%  The  in-track  formation  has  satellites  that  occupy  the  same  ground  track. 

%  The  deputy  satellite  starts  off  1  Km  behind  the  control  and  in  a 
%  slightly  different  track. 

%  This  is  accomplished  by  offsetting  their  RAAN  by 

Period_initial  =  sqrt((4*(piA2)*(Sat_SemiMajor*1000)A3)  /  (G*MassEarth)); 
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intrack_RAAN  =  (2*pi-PSat_ArgLat)/(2*pi)  *  Period_initial  *  Earth_rate; 

%  Degree  to  Rad  Conversion  and  reverse 

degtorad  =  pi/180; 
radtodeg=  180/pi; 

%  Convert  Initial  Orbit  Parameters  to 
%  Cartesian  Position  and  Velocity 
%  Uses  "posvel"  function  given  at  end  of  this  code 

[rO,vO]  =  posvel(Sat_SemiMajor,  Sat_Ecc,  Sat_Inc,  Sat_RAAN,  Sat_ArgPeri, 
Sat_ArgLat,  mu); 

[rOplane,  vOplane]  =  posvel(Sat_SemiMajor,  Sat_Ecc,  Sat_Inc,  Sat_RAAN,  Sat_ArgPeri, 
PSat_ArgLat,  mu); 

[rOtrack,  vOtrack]  =  posvel(Sat_SemiMajor,  Sat_Ecc,  Sat_Inc,  intrack_RAAN, 
Sat_ArgPeri,  PSat_ArgLat,  mu); 


stateO=[rO;vO]; 

stateplaneO=[rOplane;vOplane]; 
statetrackO=[rOtrack;vOtrack] ; 

%  Propagate  satellite  position  using  two-body  orbit  assumptions 
%  Uses  ode45  numerical  integrator  and 
%  "twobody"  function  given  at  end  of  this  code 

%  Set  max  stepsize  for  integration.  I  chose  a  number  that 
%  gives  quick  results  without  losing  too  much  accuracy 
%  over  a  10  hour  period  of  propagation 

options  =  odeset('MaxStep',50); 

[t,state]=ode45(@twobody,[t_vec],state0, options); 

staterad  =  state(:,l:3); 
statevel  =  state(:,4:6); 

statevellength  =  length(statevel); 

%Initialize  period  vector,  because  it  changes  slightly  as  drag  effects 
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%satellite 

Pvector=zeros(statevellength,  1 ); 

%Initialize  RTN  velocity  vector 
RTNstatevel=zeros(statevellength,3); 

%Initialize  true  anomaly  vector  (ta) 
ta=zeros(statevellength,  1 ); 

for  i=  1  :statcvellength 

velmag(i)=norm(statevel(i,:)); 
radmag(i)=norm(  staterad(i, : )) ; 

Pvector(i,:)=sqrt((4*(piA2)*(radmag(i)*  1000)A3)  /  (G*MassEarth)); 

%  Determine  the  true  anomaly  by  dividing  time  by  period 
%  Then  divide  remainder  by  period  and  multiply  by  2pi 
ta(i)=(rem(t_vec(i),Pvector(i))/Pvector(i))*2*pi; 

ECItoRTN  =  [(cos(Sat_ArgPeri  +  ta(i))*cos(Sat_RAAN))-(sin(Sat_ArgPeri  + 
ta(i))*cos(Sat_Inc)*sin(Sat_RAAN)),... 

cos(Sat_ArgPeri  +  ta(i))  *  sin(Sat_RAAN)  +  sin(Sat_ArgPeri  +  ta(i))  *  cos(Satlnc)* 
cos(Sat_RAAN), . . . 

sin(Sat_ArgPeri  +  ta(i))  *  sin(Sat_Inc);... 

-sin(Sat_ArgPeri  +  ta(i))  *  cos(Sat_RAAN)  -  cos(Sat_ArgPeri  +  ta(i))  *  cos(Sat_Inc)* 
sin(Sat_RAAN),. . . 

-sin(Sat_ArgPeri  +  ta(i))  *  sin(Sat_RAAN)  +  cos(Sat_ArgPeri  +  ta(i))  *  cos(Sat_Inc)  * 
cos(Sat_RAAN), . . . 

cos(Sat_ArgPeri  +  ta(i))  *  sin(Sat_Inc);... 
sin(Sat_Inc)  *  sin(Sat_RAAN),... 

-sin(Sat_Inc)  *  cos(Sat_RAAN),... 
cos(Sat_Inc)]; 

RTN  statevel(i,  :)=  ECItoRTN*  statevel(i, : 

RTNvelmag(i)=norm(RTNstatevel(i,:)); 

RTNdiff(i)=velmag(i)-RTNvebnag(i); 


end 


timehrs=t/3600; 
figure(  1 ) 

plot(timehrs,  radmag) 

title('Orbit  Radius  of  Sat#l  vs.  Time(hrs)','fontsize',16,'fontweighf,'bold') 
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figure(2) 

plot(timehrs,  velmag) 

title('Magnitude  of  Velocity  of  Sat#lvs.  Time(hrs)','fontsize',16,'fontweighf,'bold') 


%  Start  Code  for  second  satellite  with  different  drag  characteristics 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%  The  in  plane  formation 

[t2,state2]=ode45(@twobody2,[t_vec],stateplane0,options); 

staterad2  =  state2(:,l:3); 
statevel2  =  state2(:,4:6); 

statevellength2  =  length(statevel2); 


%Initialize  period  vector 
Pvector2=zeros(statevellength,  1 ); 

%Initialize  RTN  velocity  vector 
RTN  statevel2=zeros(statevellength,3 ); 

%Initialize  true  anomaly  vector  (ta) 
ta2=zeros(statevellength,  1 ); 

%Initialize  RTN  position  vector 
relpo  sRTN=zeros(statevellength,3 ); 

%Initialize  relative  position  vector 
relative_position=zeros(statevellength,3); 

for  i=l:statevellength2 

velmag2(i)=norm(statevel2(i,:)); 
radmag2  ( i)=norm(staterad2  (i, : )) ; 

Pvector2(i,:)=sqrt((4*(piA2)*(radmag2(i)*  1000)A3)  /  (G*MassEarth)); 

%  Determine  the  true  anomaly  by  dividing  time  by  period 
%  Then  divide  remainder  by  period  and  multiply  by  2pi 

ta2(i)=(rem(t_vec(i),Pvector2(i))/Pvector2(i))*2*pi; 


89 


ECItoRTN2  =  [(cos(Sat_ArgPeri  +  ta2(i))*cos(Sat_RAAN))-(sin(Sat_ArgPeri  + 
ta2(i))*cos(Sat_Inc)*sin(Sat_RAAN)),... 

cos(Sat_ArgPeri  +  ta2(i))  *  sin(Sat_RAAN)  +  sin(Sat_ArgPeri  +  ta2(i))  * 
cos(Sat_Inc)*  cos(Sat_RAAN),... 

sin(Sat_ArgPeri  +  ta2(i))  *  sin(Sat_Inc);... 

-sin(Sat_ArgPeri  +  ta2(i))  *  cos(Sat_RAAN)  -  cos(Sat_ArgPeri  +  ta2(i))  * 
cos(Sat_Inc)*  sin(Sat_RAAN),... 

-sin(Sat_ArgPeri  +  ta2(i))  *  sin(Sat_RAAN)  +  cos(Sat_ArgPeri  +  ta2(i))  * 
cos(Sat_Inc)  *  cos(Sat_RAAN),... 

cos(Sat_ArgPeri  +  ta2(i))  *  sin(Sat_Inc);... 
sin(Sat_Inc)  *  sin(Sat_RAAN),... 

-sin(Sat_Inc)  *  cos(Sat_RAAN),... 
cos(Sat_Inc)]; 

RTNstatevel2(i,:)=  ECItoRTN2*statevel2(i,:)'; 

RTNvelmag2  (i)=norm(RTN  statevel2  (i, : )) ; 

RTNveldiff2(i)=velmag2(i)-RTNvelmag2(i);  %  To  keep  the  origin  around  sat  1  do  it 
this  way 

relative_position(i,:)=staterad2(i,:)-staterad(i,:); 

relposRTN(i,:)=ECItoRTN2*relative_position(i,:)'; 

end 


timehrs2=t2/3600; 


figure(3) 

plot(timehrs2,  relposRTN(:,l)) 

title('Radial  Direction  In-Plane  'fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
figure(4) 

plot(timehrs2,  relposRTN(:,2)) 

title('Tangential  Direction  In-Plane','fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
figure(5) 

relative_position=state2(:,  1 :3)-state(:,  1 :3); 

plot3  (relative_position(: ,  1  ),relative_position( :  ,2  ),relative_position( :  ,3 )) 
xlabel('x',  'fontsize',  1 6,'fontweight','bold') 
ylabel('Y',  'fontsize',  1 6,'fontweight','bold') 
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zlabel('Z',  'fontsize',  1 6,'fontweight','bold') 
axis('equal') 

title('Relative  position’, 'fontsize',  16,'fontweight','bold') 


figure(6) 

plot(timehrs2,  relposRTN(:,3)) 

title(’Normal  Direction  In-Plane  'fontsize',  1 6,Tontwcight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 

for  i=l:statevellength2 

absdist(i)  =  norm(relposRTN(i,:)); 

end 

%  Plots  the  separation  between  the  two  in-plane  satellites 
figure(19) 

plot(timehrs,  absdist) 

title('Separation  Between  Satellites(In-Plane)',  'fontsize',  1 6, 'fontweight', 'bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 

%  Start  Code  for  satellite  in  an  in-track  orbit 
%  with  different  drag  characteristics 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%% 


[t3  ,state3  ]=ode45(@twobody3 ,  [t_vec] ,  statetrackO,  options) ; 

staterad3  =  state3(:,l:3); 
stateveB  =  state3(:,4:6); 

statevellength  =  length(stateveB); 

%Initialize  momentum  vector 
Hvector3=zeros(statevellength,3); 

%Initialize  period  vector 
Pvector3=zeros(statevellength,  1 ); 

%Initialize  RTN  velocity  vector 
RTNstatcvcl3=zcros(statcvcl  length, 3); 

%Initialize  true  anomaly  vector  (ta) 
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ta3=zeros(statevellength,  1 ); 


%Initialize  RTN  position  vector 
relposRTN3=zeros(statevellength,3); 

%Initialize  relative  position  vector 
relative_position3=zeros(statevellength,3); 

for  i=l:statevellength 
velmag3  (i)=norm(statevel3  (i, :)); 
radmag3  (i)=norm(staterad3  (i, :)) ; 

Hvector3(i,:)=cross(staterad3(i,:),  statevel3(i,:)); 

eccen3(i,:)=(l/mu)*(cross(statevel3(i,:),  Hvector3(i,:))-mu*staterad3(i,:)/radmag3(i)); 

Pvector3(i,:)=sqrt((4*(piA2)*(radmag3(i)*  1000)A3)  /  (G*MassEarth)); 

%  Determine  the  true  anomaly  by  dividing  time  by  period 
%  Then  divide  remainder  by  period  and  multiply  by  2pi 

ta3(i)=(rem(t_vec(i),Pvector3(i))/Pvector3(i))*2*pi; 

ECItoRTN3  =  [(cos(Sat_ArgPeri  +  ta3(i))*cos(intrack_RAAN))-(sin(Sat_ArgPeri  + 
ta3(i))*cos(Sat_Inc)*sin(intrack_RAAN)),... 

cos(Sat_ArgPeri  +  ta3(i))  *  sin(intrack_RAAN)  +  sin(Sat_ArgPeri  +  ta3(i))  * 
cos(Sat_Inc)*  cos(intrack_RAAN),... 

sin(Sat_ArgPeri  +  ta3(i))  *  sin(Sat_Inc);... 

-sin(Sat_ArgPeri  +  ta3(i))  *  cos(intrack_RAAN)  -  cos(Sat_ArgPeri  +  ta3(i))  * 
c o s(  S  at_Inc )  *  sin(  intr ack_RAAN ) , . . . 

-sin(Sat_ArgPeri  +  ta3(i))  *  sin(intrack_RAAN)  +  cos(Sat_ArgPeri  +  ta3(i))  * 
cos(Sat_Inc)  *  cos(intrack_RAAN),... 

cos(Sat_ArgPeri  +  ta3(i))  *  sin(Sat_Inc);... 
sin(Sat_Inc)  *  sin(intrack_RAAN),... 

-sin(Sat_Inc)  *  cos(intrack_RAAN),... 
cos(Sat_Inc)]; 

RTNstatevel3(i,:)=  ECItoRTN3*statevel3(i,:)'; 

RTNvelmag3  (i)=norm(RTN  stateveB  (i, :)); 

RTNveldiff3(i)=vebnag3(i)-RTNvelmag3(i);  %  To  keep  the  origin  around  sat  1  do  it 
this  way 

relative_position3(i,:)=staterad3(i,:)-staterad(i,:); 

relposRTN3(i,:)=ECItoRTN3*relative_position3(i,:)'; 

end 

timehrs3=t3/3600; 
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figure(7) 

plot3(relposRTN3(:,l),relposRTN3(:,2),relposRTN3(:,3)) 
title('In-Track  3-D’,  'fontsize',16,'fontweight','bold') 
xlabel('Radial  Direction  (Kilometers)') 
ylabel(’Tangential  Direction  (Kilometers)') 
zlabel('Normal  Direction  (Kilometers)') 

figure(8) 

plot(timehrs3,  relposRTN3(:,l)) 

title('Radial  Direction  In-Track',  'fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
figure(9) 

plot(timehrs3,  relposRTN3(:,2)) 

title('Tangential  Direction  In-Track',  'fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
figure(lO) 

plot(timehrs3,  relposRTN3(:,3)) 
axis([0  10-0.1  0.1]) 

title('Normal  Direction  In-Track',  'fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

%gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%% 

%  Start  code  for  satellite  in  a  circular  formation  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%% 


%  Start  determining  characteristics  of  a  satellite  in  a  circular  formation 
%  with  the  reference  satellite.  All  values  are  originally  in  an  RTN 
%  coordinate  system  with  the  reference  satellite  at  the  origin 

%  distance  (km)  between  reference  satellite  (master)  and  slave  satellite 

rOcirc  =  1 ; 

%  Mean  motion  is  needed  for  the  equations  and  we  know  that  the  two  orbits 
%  must  have  the  same  semi-major  axis  or  their  orbits  would  diverge  very 
%  quickly 
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ncirc  =  sqrt(mu/Sat_SemiMajorA3); 


%  use  equations  24-28  in  Sabol  paper  to  determine  the  other  initial 
%  conditions 

ynot  =  0; 

xnot  =  0.5; 

znot  =  xnot  *  sqrt(3); 

y  dote  ire  =  -0.0005850427710543116; 

xdotcirc  =  ynot*ncirc/2; 

zdotcirc  =  sqrt(3)*xdotcirc;  %  Use  positive  because  slave  is  at  apogee 

%  Change  the  difference  in  position  and  velocity  of  the  two  satellites  from 
%  RTN  coords  to  ECI  coords 

diffcircposORTN  =  [xnot,  ynot,  znot] 
diffcircvelORTN  =  [xdotcirc,  ydotcirc,  zdotcirc] 

%  RTNtoECI  =  ECItoRTN'; 

RTNtoECI=[l  0  0;0  0  -1;0  1  0]; 

diffcircposOECI  =  RTNtoECI  *  diffcircposORTN'; 
diffcircvelOECI  =  RTNtoECI  *  diffcircvelORTN'; 
statecircposOECI  =  diffcircposOECI  +  rO; 

statecircvelOECI  =  diffcircvelOECI  +  vO; 

[t4,state4]=ode45(@twobody4,[t_vec],statecirc0, options); 

staterad4  =  state4(:,l:3); 
statevel4  =  state4(:,4:6); 

statevellength  =  length(stateveM); 

%Initialize  period  vector 
Pvector4=zeros(statevellength,  1 ); 

%Initialize  RTN  velocity  vector 
RTNstatevel4=zeros(statevellength,3); 

%Initialize  true  anomaly  vector  (ta) 
ta4=zeros(statevellength,  1 ); 

%Initialize  RTN  position  vector 
relposRTN4=zeros(statevellength,3); 
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%Initialize  relative  position  vector 
relative_po  sition4=zeros(statevellength,  3 ) ; 


for  i=l:statevellength 

velmag4(i)=norm(statevel4(i,:)); 
radmag4(i)=norm(staterad4(i, :)) ; 

Pvector4(i,:)=sqrt((4*(piA2)*(radmag4(i)*  1000)A3)  /  (G*MassEarth)); 

%  Determine  the  true  anomaly  by  dividing  time  by  period 
%  Then  divide  remainder  by  period  and  multiply  by  2pi 

ta4(i)=(rem(t_vec(i),Pvector4(i))/Pvector4(i))*2*pi; 

crhat  =  staterad4(i,:)/norm(staterad4(i,:)); 
e_n  =  cross(staterad4(i,:),statevel4(i,:)); 
e_n_hat  =  e_n/nonn(e_n); 
e_t_hat  =  cross(e_n_hat,e_r_hat); 

ECItoRTN4  =  [e  r  hat;  e  t  hat;  e  n  hat]; 

RTNstatevel4(i,:)  =  ECItoRTN4*statevel4(i,:)'; 

relative_position4(i,:)=staterad4(i,:)-staterad(i,:); 

relposRTN4(i,:)=ECItoRTN4*relative_position4(i,:)'; 

energy4(i)  =  0.5*nonn(statevel4(i,:))A2  -  mu/nonn(staterad4(i,:)); 
end 


timehrs4=t4/3 600 ; 
figure(l  1) 

plot3(relposRTN4(:,l),relposRTN4(:,2),relposRTN4(:,3)) 
title('Circular  3-D',  'fontsize',16,'fontweighf,'bold') 
xlabel(’Radial  Direction  (Kilometers)’) 
ylabel(’Tangential  Direction  (Kilometers)’) 
zlabel(’Normal  Direction  (Kilometers)’) 

figure(12) 

plot(timehrs4,  relposRTN4(:,l)) 
axis([0  100  -0.8  0.8]) 

title(’Radial  Direction  Circular’,  ’fontsize',16,'fontweight',’bold’) 

xlabel('Time  (hours)') 

ylabel('Kilometers') 
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gtext( {'Deputy  area  =  1.5mA2’,  'Control  area  =  2mA2'}) 
figure(13) 

plot(timehrs4,  relposRTN4(:,2)) 
axis([0  100-1.8  1.8]) 

title('Tangential  Direction  Circular',  'fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
figure(14) 

plot(timehrs4,  relposRTN4(:,3)) 
axis([0  100-1.4  1.4]) 

title('Normal  Direction  Circular',  'fontsize',16,'fontweight','bold') 

xlabel('Time  (hours)') 

ylabel('Kilometers') 

gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 


end  %  main  program 

function  statedot=twobody(t, state) 

%  Equations  of  motion 
%  for  two-body  orbital  motion 
%  does  not  account  for  any  perturbations 

% 

%  States: 

%  Cartesian  Position  (states  1-3) 

%  Cartesian  Velocity  (states  4-6) 

EarthRadius  =  6378.1; 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 
mu  =  3.98601e5; 

%  Magnitude  of  the  radius  vector 
r  =  sqrt(state(  1  )A2  +  state(2)A2  +  state(3)A2); 

%  Magnitude  of  the  velocity  vector 
v  =  sqrt(state(4)A2  +  state(5)A2  +  state(6)A2); 


%  All  stuff  I'm  adding  for  drag 
%  for  a  height  of  250  km,  (Table  on  537  Vallado) 
hellp  =  r  -  Earth  Radius;  %  in  km 

hO  =  200;  %  in  km 


96 


H  =  37.105;  %inkm 

rhoO  =  2.789e-l  1;  %  in  kg/mA3 

%  compute  the  atmospheric  density  (rho)  in  kg/mA3 
%  rho  =  rhoO  *  exp(-(hellp-hO)/H); 
rho  =  rho0  *  exp(-(r-6378.1-hO)/H); 

A  =  2;  %  area  in  mA2  of  a  2m  x  2m  satellite 

Cd  =  2.2  ;  %  Nonnal  value  for  coefficient  of  drag 

mass  —  100;  %  in  kg 

B  =  (Cd*A)/mass; 


%  Equations  of  motion 

%  Since  the  satellite  is  in  a  polar  orbit  drag  will  only  be  affecting  it  in 
%  two  dimensions 

%  multiplying  the  total  Accel  of  drag  times  the  unit  vector  in  the 
%  desired  direction  gives  the  desired  Ad  component 

statedot(l:3)  =  state(4:6); 

statedot(4)  =  -(mu/rA3)*state(l)  -  (l/2*B*rho*vA2*state(4)/v); 
statedot(5)  =  -(mu/rA3)*state(2)  -  (l/2*B*rho*vA2*state(5)/v); 
statedot(6)  =  -(mu/rA3)*statc(3)  -  (l/2*B*rho*vA2*state(6)/v); 

statedot=statedot' ; 


end  %  twobody  function 

function  statedot2=twobo  dy2  (t2 ,  state2 ) 

%  Equations  of  motion 
%  for  two-body  orbital  motion 
%  does  not  account  for  any  perturbations 

% 

%  States: 

%  Cartesian  Position  (states  1-3) 

%  Cartesian  Velocity  (states  4-6) 

EarthRadius  =  6378.1; 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 
mu  =  3.98601e5; 

%  Magnitude  of  the  radius  vector 
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r2  =  sqrt(state2(l)A2  +  state2(2)A2  +  state2(3)A2); 

%  Magnitude  of  the  velocity  vector 

v2  =  sqrt(state2(4)A2  +  state2(5)A2  +  state2(6)A2); 


%  All  stuff  I'm  adding  for  drag 

%  for  a  height  of  250  km,  (Table  on  537  Vallado) 

hellp  =  r2  -  EarthRadius;  %  in  km 

hO  =  200;  %  in  km 

H  =  37.105;  %inkm 

rhoO  =  2.789e-l  1;  %  in  kg/mA3 

%  compute  the  atmospheric  density  (rho)  in  kg/mA3 
%  rho  =  rhoO  *  exp(-(hellp-hO)/H); 
rho2  =  rho0  *  exp(-(r2-6378.1-hO)/H); 

A2  =  1 .5;  %  area  in  mA2  of  a  2m  x  2m  satellite 

Cd2  =  2.2;  %  Nonnal  value  for  coefficient  of  drag 

mass2  =  100;  %  in  kg 

B2  =  (Cd2*A2)/mass2; 


%  Equations  of  motion 

%  Since  the  satellite  is  in  a  polar  orbit  drag  will  only  be  affecting  it  in 
%  two  dimensions 

%  multiplying  the  total  Accel  of  drag  times  the  unit  vector  in  the 
%  desired  direction  gives  the  desired  Ad  component 

statedot2(l:3)  =  state2(4:6); 

statedot2(4)  =  -(mu/r2A3)*state2(l)  -  (1  /2 * B 2 * rh o 2 * v 2 A 2 * s tate2 (4 )/v 2 ) ; 
statedot2(5)  =  -(mu/r2A3)*state2(2)  -  (1  /2 *  B 2 * r h o 2 * v2 A 2 * s tate2 ( 5  )/v 2 ) ; 
statedot2(6)  =  -(mu/r2A3)*state2(3)  -  (l/2*B2*rho2*v2A2*state2(6)/v2); 

statedot2=statedot2’ ; 


end  %  twobody2  function 

function  statedot3 =twobo  dy 3  (t3 ,  state3 ) 
%  Equations  of  motion 
%  for  two-body  orbital  motion 
%  does  not  account  for  any  perturbations 

% 
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%  States: 

%  Cartesian  Position  (states  1-3) 
%  Cartesian  Velocity  (states  4-6) 


EarthRadius  =  6378.1; 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 
mu  =  3.98601e5; 

%  Magnitude  of  the  radius  vector 

r3  =  sqrt(state3(l)A2  +  state3(2)A2  +  state3(3)A2); 

%  Magnitude  of  the  velocity  vector 

v3  =  sqrt(state3(4)A2  +  state3(5)A2  +  state3(6)A2); 


%  All  stuff  I'm  adding  for  drag 

%  for  a  height  of  250  km,  (Table  on  537  Vallado) 

hellp  =  r3  -  Earth  Radius;  %  in  km 

hO  =  200;  %  in  km 

H  =  37.105;  %inkm 

rhoO  =  2.789e-l  1;  %  in  kg/mA3 

%  compute  the  atmospheric  density  (rho)  in  kg/mA3 
%  rho  =  rhoO  *  exp(-(hellp-hO)/H); 
rho3  =rho0  *  exp(-(r3-6378.1-hO)/H); 

A3  =  1 .5;  %  area  in  mA2  of  a  2m  x  2m  satellite 

Cd3  =  2.2  ;  %  Nonnal  value  for  coefficient  of  drag 

mass3  =  100;  %  in  kg 

B3  =  (Cd3*A3)/mass3; 


%  Equations  of  motion 

%  Since  the  satellite  is  in  a  polar  orbit  drag  will  only  be  affecting  it  in 
%  two  dimensions 

%  multiplying  the  total  Accel  of  drag  times  the  unit  vector  in  the 
%  desired  direction  gives  the  desired  Ad  component 

statedot3(l:3)  =  state3(4:6); 

statedot3(4)  =  -(mu/r3A3)*state3(l)  -  (l/2*B3*rho3*v3A2*state3(4)/v3); 
statedot3(5)  =  -(mu/r3A3)*state3(2)  -  (l/2*B3*rho3*v3A2*state3(5)/v3); 
statedot3(6)  =  -(mu/r3A3)*state3(3)  -  (!/2*B3*rho3*v3A2*state3(6)/v3); 
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statedot3 =statedot3 ' ; 


end  %  twobody3  function 


function  statedot4=twobo  dy4(t4 ,  state4) 
%  Equations  of  motion 
%  for  two-body  orbital  motion 
%  does  not  account  for  any  perturbations 

% 

%  States: 

% 

%  Cartesian  Position  (states  1-3) 

%  Cartesian  Velocity  (states  4-6) 

% 


EarthRadius  =  6378.1; 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 
mu  =  3.98601e5; 

%  Magnitude  of  the  radius  vector 

r4  =  sqrt(state4(l)A2  +  state4(2)A2  +  state4(3)A2); 

%  Magnitude  of  the  velocity  vector 

v4  =  sqrt(state4(4)A2  +  state4(5)A2  +  state4(6)A2); 


%  All  stuff  I'm  adding  for  drag 

%  for  a  height  of  250  km,  (Table  on  537  Vallado) 

hellp  =  r4  -  Earth  Radius;  %  in  km 

hO  =  200;  %  in  km 

H  =  37.105;  %inkm 

rhoO  =  2.789e-l  1;  %  in  kg/mA3 

%  compute  the  atmospheric  density  (rho)  in  kg/mA3 
%  rho  =  rhoO  *  exp(-(hellp-hO)/H); 
rho4  =rho0  *  exp(-(r4-6378.1-hO)/H); 

A4  =1.5;  %  area  in  mA2  of  a  2m  x  2m  satellite 

Cd4  =  2.2;  %  Nonnal  value  for  coefficient  of  drag 

mass4  =  100;  %  in  kg 

B4  =  (Cd4*A4)/mass4; 
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%  Equations  of  motion 

%  Since  the  satellite  is  in  a  polar  orbit  drag  will  only  be  affecting  it  in 
%  two  dimensions 

%  multiplying  the  total  Accel  of  drag  times  the  unit  vector  in  the 
%  desired  direction  gives  the  desired  Ad  component 

statedot4(l:3)  =  state4(4:6); 

statedot4(4)  =  -(mu/r4A3)*state4(l)  -  (l/2*B4*rho4*v4A2*state4(4)/v4); 
statedot4(5)  =  -(mu/r4A3)*state4(2)  -  (l/2*B4*rho4*v4A2*state4(5)/v4); 
statedot4(6)  =  -(mu/r4A3)*state4(3)  -  ( 1  /2  *  B4  *rho4  * v4 A 2 * s tatc4 (6 )/v4 ) ; 

statedot4=statedot4'; 


end  %  twobody4  function 


function  [r,v]=posvel(a,e,inc,Omg,w,nu,mu) 

%compute  position  and  velocity  vectors  from  classical  orbital  elements 
%  [r,v]=posvel(a,  e,  inc,  Omg,  w,  nu,mu) 

%  Input  all  angles  in  radians 
%  r  and  v  output  in  IJK  unit  vectors 

% 

%  if  mu  is  not  present,  assume  canonical  units 

% 

%  CIRCULAR  ORBIT 

%  if  e=0  and  inc~=0,  w  is  not  used  and  nu  is  assumed  to  be  argument  of 
%  latitude  (uO) 

% 

%  EQUATORIAL  ORBIT 

%  if  inc=0,  Omg  is  not  used  and  w  is  assumed  to  be  Longitude  of  Periapsis 
%  (PI) 

% 

%  CIRCULAR,  EQUATORIAL  ORBIT 

%  if  e=0  and  inc=0,  Omg  and  w  are  not  used  and  nu  is  assumed  to  be  True 
%  longitude  (10) 

% 

if  nargin<6 

display('Not  enough  input  parameters') 
end 

if  nargin==6 
mu=l; 

display('Assuming  Canonical  Units  with  Mu=1.0’) 
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end 


if  e==0  &  inc~=0 
w=0.0; 

disp('Circular  Orbit,  Argument  of  Latitude  substituted  for  True  Anomaly') 
end 

if  e~=0  &  inc==0 
Omg=0.0; 

disp('Equatorial  Orbit,  Longitude  of  Periapsis  substituted  for  Argument  of  Perigee') 
end 

if  e==0  &  inc==0 
w=0; 

Omg=0; 

disp('Circular,  Equatorial  Orbit,  True  Longitude  substituted  for  True  anomaly') 
end 

%  Compute  orbit  parameter  (p) 


p=a*(l-eA2); 


%  Compute  Position  and  Velocity  Vectors  in  Perifocal  Coordinate  System 
rp  =  zeros(3,l); 

rp(l)  =  p*cos(nu)/(l+e*cos(nu)); 
rp(2)  =  p*sin(nu)/(l+e*cos(nu)); 

vp  =  zeros(3,l); 

vp(l)  =  -sqrt(mu/p)*sin(nu); 

vp(2)  =  sqrt(mu/p)*(e  +  cos(nu)); 

%  Rotate  from  Perifocal  Coordinate  System  to  IJK  system 

%  Rotation  matrix  for  -w  angle  rotation  about  axis  nonnal  to  orbit  plane 

Rw  =  [cos(-w)  sin(-w)  0;  -sin(-w)  cos(-w)  0;  0  0  1]; 

%  Rotation  Matrix  for  -inc  angle  rotation  about  first  axis  (periapsis 
%  direction  in  PQW,  which  should  be  aligned  with  node  direction  after 
%  rotation  above) 

Ri  =  [1  0  0;  0  cos(-inc)  sin(-inc);  0  -sin(-inc)  cos(-inc)]; 

%  Rotation  Matrix  for  -Omg  angle  rotation  about  third  axis  (normal  to  orbit 
%  plane,  which  should  be  aligned  with  K  after  both  rotations  above) 
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RO  =  [cos(-Omg)  sin(-Omg)  0;  -sin(-Omg)  cos(Omg)  0;  0  0  1]; 

%  Compute  Position  and  Velocity  in  IJK  by  rotation  vectors  in  PQW  coord 
r  =  RO*Ri*Rw*rp; 
v  =  RO*Ri*Rw*vp; 


end  %  posvel  function 
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Appendix  B  -  Controller  Code 


The  following  code  was  written  and  executed  using  MATLAB  version  7.0. 
function  []=orbit_propagator() 

%  Simulation  of  a  satellite  in  orbit  with  two-body  motion  and  drag 

global  variable 

variable,  time  =  []; 
variable. area  =  []; 

%  Initial  and  Final  Time  (seconds) 
tO  =  0; 

tf  =  1800000;  %  10  hours  =  36000;  1000  hours  =  3600000 
t_vec=t0:10:tf; 

%  Constants 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 
mu  =  3.98601e5; 

%  Initial  Satellite  Orbital  Parameters 
%  Orbit  Altitude  (Kilometers) 

Sat_Alt  =  250; 

EarthRadius  =  6378.1; 

%  Rotational  rate  of  Earth  (rad/s) 

Earth_rate  =  7.2722e-5;  %  corresponds  to  15  degrees/hour  ..  expressed  in  rad/s 
Sat_SemiMajor  =  Earth_Radius  +  Sat_Alt;  %  Compute  Semimajor  Axis 
%  Orbit  Eccentricity 
Sat_Ecc  =  0;  %  Eccentricity 
%  Orbit  Inclination  (radians) 
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Sat_Inc  =  pi/2; 

%  Right  Acsension  of  Ascending  Node  (radians) 

%  Angle  between  vernal  equinox  direction  and  point  at  which  the  satellite  passes 
%  the  equator  going  north 

Sat_RAAN  =  0; 

%  Argument  of  Perigee  (radians) 

%  Argument  of  Perigee  undefined  for  circular  orbit,  assume  zero  as  placeholder 
Sat_ArgPeri  =  0; 

%  Argument  of  Latitude  at  epoch  (radians) 

%  Angle  between  the  ascending  node  direction  and  the  satellite  position  vector 
%  at  the  start  time  of  the  sim 

Sat_ArgLat  =  0; 

%  Argument  of  Latitude  at  epoch  (radians) 

%  For  an  in-plane  satellite  that  is  approximately  1  km  behind  the 
%  control  satellite.  Calculation  2pi*radius  of  orbit  =  circum  of  orbit 
%  circum  of  orbit/2pi,  then  find  how  many  radians  for  1  km  separation. 

PSat_ArgLat  =  2*pi  -  0.000150873; 

%  Earth  Rotational  Position  (radians) 

%  Angle  between  prime  meridian  and  vernal  equinox  direction  at  start  time  of  sim. 

%  This  is  an  arbitrary  assignment  unless  we  want  sim  time  to  have  some  real  meaning 
%  with  respect  to  a  true  astronomical  time  system 

EarthPos  =  0; 

%  The  in-track  formation  has  satellites  that  occupy  the  same  ground  track. 

%  The  deputy  satellite  starts  off  1  Km  behind  the  control  and  in  a 
%  slightly  different  track. 

%  This  is  accomplished  by  offsetting  their  RAAN  by: 


Period_initial  =  sqrt((4*(piA2)*(Sat_SemiMajor)A3)  /  mu); 
intrack_RAAN  =  (2*pi-PSat_ArgLat)/(2*pi)  *  Period_initial  *  Earth_rate 
%  Degree  to  Rad  Conversion  and  reverse 
degtorad  =  pi/180; 
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radtodeg=  180/pi; 

%  Convert  Initial  Orbit  Parameters  to 
%  Cartesian  Position  and  Velocity 
%  Uses  "posvel"  function  given  at  end  of  this  code 

[rO,vO]  =  posvel(Sat_SemiMajor,  Sat_Ecc,  Sat_Inc,  Sat_RAAN,  Sat_ArgPeri, 
Sat_ArgLat,  mu); 

[rOplane,  vOplane]  =  posvel(Sat_SemiMajor,  Sat_Ecc,  Sat_Inc,  Sat_RAAN,  Sat_ArgPeri, 
PSat_ArgLat,  mu); 

[rOtrack,  vOtrack]  =  posvel(Sat_SemiMajor,  Sat_Ecc,  Sat_Inc,  intrack_RAAN, 
Sat_ArgPeri,  PSat_ArgLat,  mu); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%% 

%  %%%%%%%  Start  code  for  satellite  in  a  circular  formation 

%%%%%%%%%%%%% 

%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%%%%% 


%  Start  determining  characteristics  of  a  satellite  in  a  circular  formation 
%  with  the  reference  satellite.  All  values  are  originally  in  an  RTN 
%  coordinate  system  with  the  reference  satellite  at  the  origin 

%  Choose  the  separation  in  the  normal  direction  to  be  0.5km 

%  Mean  motion  is  needed  for  the  equations  and  we  know  that  the  two  orbits 
%  must  have  the  same  semi-major  axis  or  their  orbits  would  diverge  very 
%  quickly 

ncirc  =  sqrt(mu/Sat_SemiMajorA3); 

%  use  equations  24-28  in  Sabol  paper  to  determine  the  other  initial 
%  conditions  as  follows 

ynot  =  0; 

xnot  =  0.5; 

znot  =  xnot  *  sqrt(3); 

ydotcirc  =  -0.00058504277105431 16;  %  for  a  1  km  separation 

xdotcirc  =  ynot*ncirc/2; 

zdotcirc  =  sqrt(3)*xdotcirc;  %  Use  positive  because  slave  is  at  apogee 
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diffcircposORTN  =  [xnot,  ynot,  znot]; 
diffcircvelORTN  =  [xdotcirc,  ydotcirc,  zdotcirc]; 

%  RTNtoECI  =  ECItoRTN'; 

RTNtoECI=[l  0  0;0  0  -1;0  1  0]; 

diffcircposOECI  =  RTNtoECI  *  diffcircposORTN'; 
diffcircvelOECI  =  RTNtoECI  *  diffcircvelORTN'; 

rOcircular  =  diffcircposOECI  +  rO; 
vOcircular  =  diffcircvelOECI  +  vO; 


%  rOcircular=[xnot;ynot;znot]; 

%  vOcircular=  [xdotcirc ;ydotcirc ;zdotcirc] ; 


%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

!!!!!!!!!!% 

%  Choose  the  type  of  formation  you  want  to  analyze  here,  and  comment  out 
%  the  othere  initial  conditions  for  the  different  formations  %%%%%%%%%% 
State0=[r0;v0;r0plane;v0plane;0];  %  Set  to  look  at  two  in-plane  satellites 
%state0=[r0;v0;r0traek;v0traek;0];  %  Set  to  look  at  two  in-track  satellites 
%state0=[r0;v0;r0eircular;v0circular;0];  %Set  to  look  at  two  circular  formation  sats 

%  Note  there  is  a  13th  state  above.  This  state  is  going  to 
%  integrate  the  difference  in  orbital  energy  and  is  part  of  the  controller 
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

%  Propagate  satellite  position  using  two-body  orbit  assumptions 
%  Uses  ode45  numerical  integrator  and 
%  "twobody_wdrag"  function  given  at  end  of  this  code 

%  Set  max  stepsize  for  integration.  I  chose  a  number  that 
%  gives  quick  results  without  losing  too  much  accuracy 
%  over  a  10  hour  period  of  propagation 

options  =  odeset('MaxStep’,50); 

[t,state]=ode45(@twobody_wdrag,[t_vec],state0, options); 

statelength  =  length(state); 

stateradl  =  state(:,l:3)’; 
statevell  =  state(:,4:6)'; 


107 


staterad2  =  state(:,7:9)'; 
statevel2  =  state(:,10:12)'; 

d_pos  =  zeros(statelength,3); 
d_vel  =  zeros(statelength,3); 

for  i=l  :statelength 

crhat  =  stateradl(:,i)/norm(stateradl(:,i)); 
e_n  =  cross(stateradl(:,i),statevell(:,i)); 
e_n_hat  =  e_n/nonn(e_n); 
cthat  =  cross(e_n_hat,e_r_hat); 

ECItoRTN  =  [e  r  hat,  e  t  hat,  e  n  hat]'; 

%  RTNstatevel(i,:)=  ECItoRTN*statevel(i,:)'; 
%  RTNvelmag(i)=norm(RTNstatevel(i,:)); 

%  RTNdiff(i)=velmag(i)-RTNvelmag(i); 


d_pos(i,:)  =  ECItoRTN*(staterad2(:,i)  -  stateradl(:,i)); 
d_vel(i,:)  =  ECItoRTN*(statevel2(:,i)  -  statevell(:,i)); 

d_energy(i)  =  0.5  *  norm(staterad2(: , i ))A2  -  mu/norm(staterad2(:,i))  -  ... 
0.5*norm(  staterad  1  ,  i  ))A2  +  mu/norm(stateradl(:,i)); 


end 

timehrs=t/3600; 
figure(  1 ) 

plot(timehrs,d_energy) 

figure(2) 

plot(d_pos(:,l),d_pos(:,2)) 

%axis([-0.001  0.001  -1.001  -.999])  %In-Plane  parameters 
title('Tangential  vs.  Radial  Separation  In-Plane  'fontsize',16,'fontweighf,'bold') 
xlabel('Radial  Separation  (Km)’) 
ylabel('Tangential  Separation  (Km)’) 

figure(3) 

plot(timehrs,d_pos) 


figure(4) 

plot3(d_pos(:,l),d_pos(:,2),d_pos(:,3)) 

xlabel('Radial  Separation  (Km)’) 
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ylabel('Tangential  Separation  (Km)') 
zlabel('Nonnal  Separation  (Km)') 


figure(5) 

plot(timehrs,  d_pos(:,l)) 

title('Radial  Direction  In-Plane  ',  'fontsize',16,'fontweight','bold') 
axis([0  500  -10e-4  8e-4])  %In-Plane  parameters 
xlabel('Time  (hours)') 
ylabel('Kilometers') 

gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
gtext({'kP  =  15','kl  =  2'}) 

figure(6) 

plot(timehrs,  d_pos(:,2)) 

title('Tangential  Direction  In-Plane  ',  'fontsize',16,'fontweight','bold') 
axis([0  500  -1.7  -0.5])  %In-Plane  parameters 
xlabel('Time  (hours)') 
ylabel('Kilometers') 

gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
gtext({'kP  =  15','kl  =  2'}) 

figure(7) 

plot(timehrs,  d_pos(:,3)) 

title('Normal  Direction  In-Plane  ',  'fontsize',16,'fontweight','bold') 

%axis([0  500 -1.5  1.5]) 
xlabel('Time  (hours)') 
ylabel('Kilometers') 

gtext( {'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
gtext({'kP  =  15','kl  =  2'}) 

figure(8) 

plot(variable. time/3600,  variable. area) 

title('Cross  Sectional  Area  vs.  Time',  'fontsize',16,'fontweight','bold') 
xlabel('Time  (hours)') 
ylabel('Square  Meters') 

gtext({ 'Initial  Conditions',  'Deputy  area  =  1.5mA2',  'Control  area  =  2mA2'}) 
gtext({'kP  =  15','kl  =  2'}) 


end  %  main  program 


function  statedot=twobody_wdrag(t, state) 
%  Equations  of  motion 
%  for  two-body  orbital  motion 
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%  and  adds  drag  plus  a  drag  controller 

% 

%  States: 

% 

%  Satellite  1  Position  (states  1-3) 

%  Satellite  1  Velocity  (states  4-6) 

%  Satellite  2  Position  (states  7-9) 

%  Satellite  2  Velocity  (states  10-12) 

%  State  13  used  for  controller 


global  variable 

%  Earth  Gravitational  Parameter  (kmA3/secA2) 
mu  =  3.98601e5; 

Earth_Radius  =  6378.1;  %  kilometers 


%  Magnitude  of  the  radius  vector  (sat  1) 
rl  =  sqrt(state(l)A2  +  state(2)A2  +  state(3)A2); 

%  Magnitude  of  the  velocity  vector  (satl) 
vl  =  sqrt(state(4)A2  +  state(5)A2  +  state(6)A2); 

%  Magnitude  of  the  radius  vector  (sat  2) 
r2  =  sqrt(state(7)A2  +  state(8)A2  +  state(9)A2); 

%  Magnitude  of  the  velocity  vector  (sat  2) 

v2  =  sqrt(state(10)A2  +  state(l  1)A2  +  state(12)A2); 

%  Drag  Modeling 

%  for  a  height  of  250  km,  (Table  on  537  Vallado) 

%  atmospheric  density  is  based  on  satellite  1  altitude 
hellp  =  rl  -  Earth  Radius;  %  in  km 
hO  =  200;  %  in  km 

H  =  37.105;  %inkm 

rhoO  =  2.789e-l  1;  %  in  kg/mA3 

%  Drag  Modeling 

%  for  a  height  of  250  km,  (Table  on  537  Vallado) 

%  atmospheric  density  is  based  on  satellite  2  altitude 
hellp2  =  r2  -  Earth_Radius;  %  in  km 
hO  =  200;  %  in  km 

H  =  37.105;  %inkm 

rhoO  =  2.789e-l  1;  %  in  kg/mA3 
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%  compute  the  atmospheric  density  (rho)  in  kg/mA3 
%  rho  =  rhoO  *  exp(-(hellp-hO)/H); 
rho  =  rhoO  *  exp(-(rl-6378.1-hO)/H); 


%  compute  the  atmospheric  density  (rho)  in  kg/mA3 
%  rho  =  rhoO  *  exp(-(hellp-hO)/H); 
rho2  =  rhoO  *  exp(-(r2-6378.1-hO)/H); 


%  Compute  Ballistic  Coefficient  for  satelite  one 


%  nominal  area  in  mA2  of  a  micro  satellite 
%  delta-area  (from  damaged  satellite,  for  example) 
%  Guesstimate  for  coefficient  of  drag 
%  in  kg 


AO  =  1.5; 
dAl  =  0.5; 

Cd  =  2.2; 
mass  =100; 

Al=A0+dAl; 

B1  =  (Cd*Al)/mass;  %  ballistic  coefficient 


% - Controller  - 

% 

%  Simple  controller  based  on  the  period  (energy)  difference  between  the  two 
%  satellites.  Computes  Specific  Mechanical  Energy  for  both  spacecraft  and 
%  then  assigns  an  Area  for  satellite  two  that  is  proportional  to  the  energy 
%  difference.  If  energy  of  satellite  two  is  greater  than  satellite  one, 

%  increase  the  Area  of  satellite  two,  increasing  drag,  and  lowering  energy. 

%  Compute  Energy  Difference 

El  =  0.5*vlA2  -  mu/rl; 

E2  =  0.5*v2A2  -  mu/r2; 

delta_E  =  E2-E1; 

%  Set  second  sat  area 

kP=15;  %1.5e5;  %  initial  value  %  Proportional  Controller  Gain - 

kl  =  2;  %  2e4;  %initial  value  %  Integral  Controller  Gain - 

A2  =  AO  +  kP*delta_E  +  kl*state(13); 

variable.time  =  [variable. time  t]; 
variable. area  =  [variable. area  A2]; 
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B2  =  (Cd*A2)/mass;  %  ballistic  coefficient 


%  Equations  of  motion 

%  Since  the  satellite  is  in  a  polar  orbit  drag  will  only  be  affecting  it  in 
%  two  dimensions 

%  multiplying  the  total  Accel  of  drag  times  the  unit  vector  in  the 
%  desired  direction  gives  the  desired  Ad  component 

statedot(l:3)  =  state(4:6); 

statedot(4)  =  -(mu/rlA3)*state(l)  -  (l/2*Bl*rho*vl*state(4)); 
statedot(5)  =  -(mu/rlA3)*state(2)  -  (l/2*Bl*rho*vl*state(5)); 
statedot(6)  =  -(mu/rlA3)*state(3)  -  (l/2*Bl*rho*vl*state(6)); 

statedot(7:9)  =  state(10:12); 

statedot(lO)  =  -(mu/r2A3)*state(7)  -  (l/2*B2*rho*v2*state(10)); 
statedot(l  1)  =  -(mu/r2A3)*state(8)  -  (l/2*B2*rho*v2*state(l  1)); 
statedot(12)  =  -(mu/r2A3)*state(9)  -  (l/2*B2*rho*v2*state(12)); 
statedot(13)  =  delta_E;  %  This  state  is  needed  for  the  integral  control 

statedot=statedot' ; 

end  %  twobody_wdrag  function 


function  [r,v]=posvel(a,e,inc,Omg,w,nu,mu) 

%compute  position  and  velocity  vectors  from  classical  orbital  elements 
%  [r,v]=posvel(a,  e,  inc,  Omg,  w,  nu,mu) 

%  Input  all  angles  in  RADIANS 
%  r  and  v  output  in  IJK  unit  vectors 

% 

%  if  mu  is  not  present,  assume  canonical  units 

% 

%  CIRCULAR  ORBIT 

%  if  e=0  and  inc~=0,  w  is  not  used  and  nu  is  assumed  to  be  argument  of 
%  latitude  (uO) 

% 

%  EQUATORIAL  ORBIT 

%  if  inc=0,  Omg  is  not  used  and  w  is  assumed  to  be  Longitude  of  Periapsis 
%  (PI) 

% 

%  CIRCULAR,  EQUATORIAL  ORBIT 

%  if  e=0  and  inc=0,  Omg  and  w  are  not  used  and  nu  is  assumed  to  be  True 
%  longitude  (10) 

% 
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if  nargin<6 

display('Not  enough  input  parameters') 
end 

if  nargin==6 
mu=l; 

display('Assuming  Canonical  Units  with  Mu=1.0') 
end 

degtorad  =  pi/180; 

if  e==0  &  inc~=0 
w=0.0; 

disp('Circular  Orbit,  Argument  of  Latitude  substituted  for  True  Anomaly') 
end 

if  e~=0  &  inc==0 
Omg=0.0; 

disp('Equatorial  Orbit,  Longitude  of  Periapsis  substituted  for  Argument  of  Perigee') 
end 

if  e==0  &  inc==0 
w=0; 

Omg=0; 

disp('Circular,  Equatorial  Orbit,  True  Longitude  substituted  for  True  anomaly') 
end 

%  Compute  orbit  parameter  (p) 


p=a*(l-eA2); 


%  Compute  Position  and  Velocity  Vectors  in  Perifocal  Coordinate  System 
rp  =  zeros(3,l); 

rp(l)  =  p*cos(nu)/(l+e*cos(nu)); 
rp(2)  =  p*sin(nu)/(l+e*cos(nu)); 

vp  =  zeros(3,l); 

vp(l)  =  -sqrt(mu/p)*sin(nu); 

vp(2)  =  sqrt(mu/p)*(e  +  cos(nu)); 

%  Rotate  from  Perifocal  Coordinate  System  to  IJK  system 
%  Rotation  matrix  for  -w  angle  rotation  about  axis  nonnal  to  orbit  plane 
Rw  =  [cos(-w)  sin(-w)  0;  -sin(-w)  cos(-w)  0;  0  0  1]; 
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%  Rotation  Matrix  for  -inc  angle  rotation  about  first  axis  (periapsis 
%  direction  in  PQW,  which  should  be  aligned  with  node  direction  after 
%  rotation  above) 


Ri  =  [1  0  0;  0  cos(-inc)  sin(-inc);  0  -sin(-inc)  cos(-inc)]; 

%  Rotation  Matrix  for  -Omg  angle  rotation  about  third  axis  (normal  to  orbit 
%  plane,  which  should  be  aligned  with  K  after  both  rotations  above) 

RO  =  [cos(-Omg)  sin(-Omg)  0;  -sin(-Omg)  cos(Omg)  0;  0  0  1]; 

%  Compute  Position  and  Velocity  in  IJK  by  rotation  vectors  in  PQW  coord 

r  =  RO*Ri*Rw*rp; 

v  =  RO*Ri*Rw*vp; 

end  %  posvel  function 


114 


Bibliography 


Adams,  C.,  A.  Robertson,  K.  Zimmerman,  and  J.P.  How,  “Technologies  for  Spacecraft 
Formation  Flying,”  Proceedings  of  ION  GPS-96  Conference.  Kansas  City  MO, 
(September  1996). 

Bauer,  Frank  H.,  Kate  Hartman,  Jonathon  P.  How,  John  Bristow,  David  Weidow,  and 
Franz  Busse.  “Enabling  Spacecraft  Formation  Flying  through  Spaceborne  GPS  and 
Enhanced  Automation  Technologies,”  ION  GPS-99  Conference,  369-383.  Nashville  TN, 
(September  1999). 

Carpenter,  J.,  Jesse  Leitner,  Richard  Burns,  and  David  Folta.  “Benchmark  Problems  for 
Spacecraft  Formation  Flying,”  Proceedings  of  the  2003  Guidance,  Navigation,  and 
Control  Conference  and  Exhibit,  Austin  TX  (August  2003). 

Carpenter,  J.R.,  David  Folta,  and  D.A  Quinn,  “Integration  of  Decentralized  Linear- 
Quadratic-Gaussian  Control  into  GSFC’s  Universal  3-D  Autonomous  Formation  Flying 
Algorithm.”  AIAA  GNC,  (August  1 999). 

Cobb,  Richard  G.  Air  Force  Institute  of  Technology  Professor,  Dayton  OH.  Personal 
Interview.  22  Nov  2005. 

Corazzini,  T.,  A.  Robertson,  J.C.  Adams,  A.  Hassibi,  and  J.P.  How.  “GPS  Sensing  for 
Spacecraft  Formation  Flying,”  Proceedings  of  the  ION-GPS-97  Conference,  (September 
1997). 

Clohessy  W.H.  and  R.S.  Wiltshire.  “Terminal  Guidance  System  for  Satellite 
Rendezvous,”  Journal  of  the  Aerospace  Sciences  27:9  (September  1960). 

de  Selding,  Peter  B.  “CNES  Readies  Essaim  Satellites  for  December  Launch,”  Space 
News  (October  2004).  14  October  2005 

http://www.space.com/spacenews/archive04/cnesarch  100404.html. 

Ferguson,  Phillip,  Franz  Busse,  Brian  Enberg,  Jonathan  How,  Michael  Tillerson,  Nick 
Pohlman,  Arthur  Richards,  and  Robert  Twiggs.  “Formation  Flying  Experiments  on  the 
Orion-Emerald  Mission,”  Proceedings  of  the  AIAA  2001  Space  Conference  and 
Exposition.  28-30.  Albuquerque  NM  (August  2001). 

Folta,  David,  John  Bristow,  Albin  Hawkins,  and  Greg  Dell.  “  NASA’s  Autonomous 
Formation  Flying  Technology  Demonstration,  Earth  Observing- 1  (EO-1),”  Proceedings 
of  the  International  Symposium  on  Formation  Flying,  Toulouse,  France  (October  2002). 

Francis,  Scott,  Samuel  Stein,  Thomas  Celano,  Jeremy  Warriner,  Jesse  Leitner,  Michael 
Moreau,  Richard  Burns,  Robert  Nelson,  and  A1  Gifford.  “Analytical  Tools  for  Clocks  in 


115 


Space.”  Proceedings  of  the  2003  IEEE  Frequency  Control  Symposium.  Tampa  FL  (May 
2003). 

Gaposhkin,  E.M.  and  A.  J.  Coster.  “Evaluation  of  New  Parameters  for  Use  in 
Atmospheric  Models,”  Paper  AAS-87-555.  Proceedings  of  the  1987  AAS/AIAA 
Astrodynamics  Conference.  San  Diego  CA:  AAS  Publication  Office,  1988 

Gill,  Eberhand.  “First  Results  From  a  Hardware-in-the-loop  Demonstration  of  Closed- 
loop  Autonomous  Formation  Flying,”  Advances  in  the  Astronautical  Sciences,  1 13: 361- 
375.  (2003). 

How,  J.P.,  R.  Twiggs,  D.  Weidow,  K.  Hartman,  and  F.  Bauer.  “Orion:  A  Low-Cost 
Demonstration  of  Formation  Flying  in  Space  Using  GPS,”  Proceedings  of  the  AIAA 
Astrodynamics  Specialists  Conference.  (August  1998). 

Kang,  W.,  Andy  Sparks,  and  Siva  Banda.  “Multi- Satellite  Formation  and 
Reconfiguration,”  Proceedings  of  the  2000  American  Control  Conference.  379-383. 
Chicago  IL,  (June  2000). 

Kitts,  Christopher,  Robert  Twiggs,  and  Jonathon  How.  “Emerald:  A  Low-Cost  Spacecraft 
Mission  for  Validating  Fonnation  Flying  Technologies,”  Proceedings  of  the  IEEE 
Aerospace  Conference.  2:  217-226  (March  1999) 

Larson,  Wiley  J.  and  James  Wertz.  Space  Mission  Analysis  and  Design.  Torrance  CA: 
Microcosm  Press,  1991. 

Leitner,  Jesse.  “Fonnation  Flying:  The  Future  of  Remote  Sensing  from  Space,” 
Proceedings  of  the  International  Symposium  on  Space  Flight  Dynamics.  Munich 
Germany  (December  2004). 

Leonard,  C.L.,  W.M.  Hollister,  and  E.V.  Bergmann.  “Orbital  Forma tionkeeping  with 
Differential  Drag,”  Journal  of  Guidance,  Control,  and  Dynamics,  12:  108-113  (January 
1989). 

MacDowall,  R.,  S.  Bale,  L.  Gopalswarmy,  D.  Jones,  M.  Kaiser,  J.  Kasper,  M.  Reiner,  and 
K.  Weiler.  “Solar  Imaging  Radio  Array  (SIRA):  A  Multi-Spacecraft  Mission,” 
Proceedings  of  SPIE.  284-292  (2005). 

Mohammed,  John  L.,  “Mission  Planning  for  a  Formation-Flying  Cluster,”  Proceedings  of 
the  14th  Annual  FLAIRS  Conference,  58-62,  Key  West  FL,  2001. 

Morring,  Frank  Jr.  “Formation  Flying,”  Aviation  Week  and  Space  Technology,  163:  38 
(September  2005). 

Ogata,  Katsuhiko.  Modern  Control  Engineering.  Englewood  Cliffs  CA:  Prentice-Hall 
Inc.  151-156.  1970. 


116 


Richelson,  Jeffrey  T.  “U.S.  Satellite  Imagery,  1960-1999,”  National  Security  Archive 
Electronic  Briefing  Book  No.  13,  (April  1999). 

Roux,  A.,  “Clusters  Regroup  for  Relaunch,”  Aerospace  America,  36:  48-51  (August 
1998). 

Sabol,  Chris,  Rich  Burns,  and  Craig  McLaughlin.  “Satellite  Formation  Flying  Design  and 
Evolution,”  Journal  of  Spacecraft  and  Rockets,  38:  270-278  (March-April  2001). 

Scharf,  Daniel  P.,  Fred  Y.  Hadaegh,  and  Scott  R.  Ploen.  “A  Survey  of  Spacecraft 
Formation  Flying  Guidance  and  Control  (Part  I),”  Proceedings  of  the  American  Control 
Conference,  2:  1733-1739  (June  2003). 

Storz,  Mark.  “High  Accuracy  Drag  Model  (HASDM),”  Proceedings  of  the  AIAA/AAS 
Astrodynamics  Specialist  conference  and  exhibit,  Monterey  CA  2002. 

Vallado,  David  A.  Fundamentals  of  Astrodynamics  and  Applications.  El  Segundo  CA: 
Microcosm  Press,  2001. 

Wiesel,  William  E.  Spaceflight  Dynamics.  Boston  MA:  Irwin/McGraw-Hill,  1997. 

Zetocha,  Paul,  “Satellite  Cluster  Command  and  Control,”  Proceedings  of  the  2000  IEEE 
Aerospace  Conference.  7:  49-54.  2000. 


117 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  074-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data 
sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other 
aspect  of  the  collection  of  information,  including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information 
Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other 
provision  of  law,  no  person  shall  be  subject  to  an  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 
PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

23-03-2006  Master’s  Thesis 


4.  TITLE  AND  SUBTITLE 

Characterizing  and  Controlling  the  Effects  of  Differential 
Drag  on  Satellite  Formations 


3.  DATES  COVERED  (From  -  To) 

August  2005  -  March  2006 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


Wedekind,  James  T.,  Captain,  USAF 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


7.  PERFORMING  ORGANIZATION  NAMES(S)  AND  ADDRESS(S) 

Air  Force  Institute  of  Technology 

Graduate  School  of  Engineering  and  Management  (AFIT/EN) 
2950  Hobson  Way,  Building  640 
WPAFB  OH  45433-8865 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFRF/VS 

Kirtland  AFB,  NM  87117 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


5f.  WORK  UNIT  NUMBER 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFIT/GSS/ENY/06-M 14 


10.  SPONSOR/MONITOR’S 
ACRONYM(S) 

AFRL/VS 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


14.  ABSTRACT 

The  ability  to  fly  satellites  in  close  formations  represents  a  capability  that  could  revolutionize  the  way  satellite 
missions  are  designed  in  the  future.  This  study  examines  three  of  the  primary  formation  flying  designs  and 
characterizes  the  affect  that  an  anomalous  satellite  with  a  slightly  different  cross-sectional  area  would  have  on  the 
stability  of  the  formation.  Following  the  characterization  of  the  effects  a  controller  is  implemented  to  mitigate  the 
cross-sectional  area  differences  between  the  satellites.  With  the  addition  of  a  straightforward  controller,  small  cross- 
sectional  area  differences  can  be  mitigated  and  corrected  such  that  the  satellites  will  remain  in  close  proximity  and  in 
some  cases  the  formation  will  remain  stable. 


15.  SUBJECT  TERMS 

Satellite  Constellations,  Aerodynamic  Drag,  Maneuvering  Satellites,  and  Orbiting  Satellites 


118 


